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ABSTRACT 


Systematic  investigation  is  made  of  effects  of  kinematic  assumptions  and 
finite  element  approximations  in  the  context  of  nonlinear  flexible 
multibody  dynamics.  Two  nonlinear  beam  finite  elements  are  consistently 
derived  from  virtual  work  principle  using  BemoullitEuler  and 
Timoshenko  beam  kinematics.  Initial  assessment  is  made  by  studying 
convergence  properties  of  element  formulations  with  eigenvalue  problems: 
free  vibration,  static  buckling,  and  dynamic  buckling.  Equations  of  motion 
are  derived  for  rigid  central  body  with  flexible  appendage  using  virtual 
work  principle.  Virtual  work  principle  allows  natural  and  consistent 
discretization  of  flexible  appendage  using  finite  element  method. 
Nonlinearities  in  flexibility  are  explored  through  dynamics  examples  using 
beam  finite  elements.  Application  of  dynamics  formulation  is  made  to  a 
realistic  scenario  involving  space  shuttle  remote  manipulator  arm  with 
attached  payload.  Contribution  of  nonlinear  theory,  in  both  formulation 
and  solution,  is  assessed. 
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Chapter  1 

Introduction 

1.1  Purpose/Obiective  of  Present  Work 

Increasingly,  formulations  for  flexible  multibody  dynamical  systems 
are  employing  the  finite  clement  method  in  the  discretization  of  the  flexible 
domain.  Embedded  in  the  finite  element  method  are  assumptions  regarding 
the  assumed  displacement  field,  and  additional  approximations  such  as  mass 
lumping  and  reduced  integration  over  the  spatial  domain. 

A  study  is  conducted  to  address  the  application  of  finite  element 
discretization  in  flexible  multibody  dynamics.  The  virtual  work  principle 
is  chosen  as  the  basis  for  derivation  of  the  equations  of  motion  for  a  simple 
class  of  satellite-type  vehicles.  The  reasons  for  choosing  the  virtual  work 
principle  are  threefold: 

•  an  integral  representation  of  the  governing  equations  of 
motion  is  embedded  in  the  virtual  work  principle, 

•  the  virtual  work  principle  allows  decomposition  of  dynamic 
system  into  rigid  and  flexible  portions, 

•  the  virtual  work  principle  is  a  basis  for  the  finite  element 
method. 

Kinematic  assumptions  and  finite  element  approximations  are 
investigated  in  a  dynamics  context  through  a  series  of  eigenvalue  problems. 
The  approach  emphasizes  understanding  of  the  behavior  of  consistently 
derived  finite  elements  rather  than  demonstrating  one  formulation  over 
another.  The  consistent  derivation  of  nonlinear  finite  elements  allows  an 
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assessment  of  such  ‘inconsistent’  assumptions  as  lumped  mass  and  reduced 
integration.  Such  in-depth  study  of  the  behavior  of  finite  elements  in 
dynamics  is  not  widely  available  in  the  literature. 

Multibody  dynamic  formulations  are  inherently  nonlinear  due  to  the 
large  rotations  of  reference  frames  in  inertial  space.  When  coupled  with 
the  possibility  of  nonlinear  flexibility,  the  importance  of  one  effect 
compared  to  another  is  unclear.  The  available  literature  is  not  clear  in  the 
meaning  of  ‘nonlinear’  solutions,  since  nonlinearity  arises  from  both 
inertial  and  flexibility  considerations.  Researchers  have  not  taken  an 
engineering  perspective;  they  have  not  made  explicit  statements  regarding 
the  most  important  effects,  i.e.  which  effects  are  essential  to  capture  the 
physics  of  the  problem.  The  nonlinear  equations  of  motion  developed  in 
this  thesis  are  applied  to  dynamic  problems  with  an  eye  toward 
understanding  the  separate  effects  of  forcing  terms  and  the  relative 
importance  of  nonlinear  flexibility. 


1.2  Thesis  Overview 

Chapter  2  derives  exact  integral  equations  of  motion  for  a  vehicle 
composed  of  central  rigid  body  with  a  rigidly  attached  flexible  appendage. 
Assumptions  of  lumped  mass  and  lumped  mass/inertia  are  employed  to 
yield  equations  of  motion  which  are  suitable  for  time  integration  using  an 
implicit  scheme.  Solution  of  these  equations  allows  assessment  of  the 
influence  of  centrifugal  and  Coriolis  forcing  terms  and  importance  of 
nonlinear  flexibility  (geometric  stiffness). 


17 


Chapter  1:  Introduction 


In  chapter  3,  two  nonlinear  beam  finite  elements  are  consistently 
derived  from  the  virtual  work  principle.  Bemoulli-Euler  and  Timoshenko 
beam  kinematics  are  employed  to  give  isoparametric  beam  finite  elements 
with  Cl  and  continuity,  respectively.  It  is  shown  that  consistent 
derivation  produces  higher  order  stress  resultants  in  the  geometric  stiffness 
matrix,  which  are  generally  ignored.  For  completeness,  consistent  nodal 
loads  are  also  derived. 

Systematic  assessment  of  beam  finite  elements  is  made  in  chapter  4 
through  eigenvalue  problems  which  include  free  vibration,  and  static  and 
dynamic  buckling.  Effects  of  beam  kinematics  and  finite  element 
assumptions  are  explored  in  detail  and  compared  with  the  analytical 
solution  from  Bemoulli-Euler  beam  theory.  Results  provide  insight  into 
modelling  considerations  typical  of  scenarios  arising  in  flexible  multibody 
dynamics. 

Nonlinear  dynamic  problems  are  addressed  in  chapter  5  using  the 
lumped  mass/inertia  formulation  of  chapter  2.  Two  examples  are 
considered:  the  beam  spin-up  problem,  and  a  realistic  application 
involving  the  orbiter/remote  manipulator  system  (RMS).  System  response 
is  evaluated  in  terms  of  the  level  of  discretization,  and  the  contribution  of 
nonlinear  effects  is  addressed.  Nonlinear  effects  due  to  rigid  body  motion 
coupled  with  flexible  deformations  are  differentiated  from  flexible 
nonlinearity  associated  with  tangent  stiffness  matrix.  Gross  overall  motion 
of  the  system  and  flexible  natural  frequency  are  considered  independently. 

Six  appendices  are  included  for  completeness.  Stress  resultants 
which  arise  in  beam  elements  are  considered  in  Appendix  A.  In  Appendix 
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B,  the  Newmark  integration  scheme  is  derived  for  linear  systems  and  an 
incremental  form  with  Newton/Raphson  iteration  is  derived  for  solution  of 
nonlinear  systems.  Introduction  to  the  explicit  Runge-Kutta  integration 
scheme  is  also  given.  Convergence  data  is  tabulated  for  future  reference  in 
Appendix  C.  Details  of  the  RMS  finite  element  model  are  given  in 
Appendix  D.  In  Appendix  E,  the  Timoshenko  frequency  equation 
(cantilevered  boundary  conditions)  is  solved,  and  the  natural  frequency 
convergence  plots  repeated  with  alternative  normalization.  Finally,  a  brief 
overview  of  Gauss  quadrature  is  given  in  Appendix  F. 


1.3  Literature  Review 

1.3.1  Beam  Theorv/Finite  Eiements 

Classical  Bemoulli-Euler  beam  theory  is  known  to  overpredict  the 
natural  frequencies  for  higher  modes  of  vibration.  It  also  tends  to 
overpredict  natural  frequencies  for  all  modes  for  thick  beams 
(length/thickness  <  10).  The  first  problem  was  alleviated  by  Rayleigh  [1], 
who  introduced  rotatory  inertia  of  the  beam  cross-sections.  An  additional 
modification  was  introduced  by  Timoshenko  [2,  3],  allowing  description  of 
cross-section  and  neutral  axis  rotation  by  independent  angles,  thus  allowing 
the  beam  to  undergo  shearing  deformation. 

The  literature  has  focused  on  Timoshenko  beam  theory,  from  both 
analytical  and  finite  element  standpoints.  The  partial  differential  equations 
resulting  from  Timoshenko’s  theory  are  difficult  to  solve  for  anything  but 
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prismatic  beams  with  simple  boundary  conditions.  Huang  [4]  derives 
frequency  and  normal  mode  equations  for  uniform,  isotropic  beams  with 
simple  boundary  conditions.  Noting  the  difficulty  of  solving  the  frequency 
equations,  he  introduces  the  ‘frequency  chart’,  which,  for  a  given  set  of 
beam  parameters,  provides  a  correction  factor  to  be  applied  to  the 
Bemoulli-Euler  solution  for  natural  frequency.  Frequency  charts  provide 
a  quick  and  convenient  measure  of  the  influence  of  rotatory  inertia  and 
shearing  deformation.  Leckie  and  Lindberg  [5]  studied  the  effect  of 
lumped  mass  assumptions  on  beam  natural  frequencies  using  finite 
difference  expressions. 

In  the  finite  element  literature,  much  emphasis  has  been  placed  on 
the  development  of  higher  order  Timoshenko  beam  elements.  Higher 
order  elements  [6,  7,  8]  were  necessary  in  order  to  satisfy  all  geometric 
and  natural  boundary  conditions  of  a  Timoshenko  beam.  The  simplest 
shear  deformable  beam  possible  was  introduced  by  Hughes  et.  al.  [9]. 
Convergence  and  accuracy  were  demonstrated  for  static  problems.  Shear 
locking  was  avoided  by  use  of  selective  reduced  integration.  Consistent 
assessment  of  finite  element  approximations  for  dynamics  has  not  been 
undertaken  in  the  literature. 

1.3.2  Multibodv  Dynamics 

Review  and  chronology  of  the  rigid  and  flexible  multibody  dynamics 
literature  is  widely  available  [10,  11,  12]  and  will  not  be  repeated  here. 
The  literature  can  be  further  partitioned  according  to  the  intended 
application.  Mello  [10]  separated  the  literature  into  the  following  groups: 
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spacecraft,  mechanisms,  and  robotics.  There  exists  a  body  of  literature 
whose  multibody  formulations  are  demonstrated  using  beam  finite 
elements.  A  selection  of  these  will  be  addressed. 

Ryan  [13]  investigated  a  deficiency  of  conventional  general  flexible 
multibody  programs  (such  as  DISCOS,  NBOD,  TREETOPS,  ALLFLEX), 
which  fail  to  correctly  capture  the  stiffening  effect  of  rapidly  spinning 
systems.  He  extended  the  assumed  modes  formulation  and  demonstrated 
the  new  theory  by  application  to  a  deployment  maneuver  and  the  beam 
spin-up  problem.  Simo  &  Vu-Quoc  [14,  15]  developed  equations  of  motion 
for  a  flexible  beam  undergoing  large  overall  motions.  In  [15],  the  spin-up 
problem  is  addressed  among  other  examples.  Quadratic  beam  finite 
elements  are  used  in  the  discretization  of  the  flexible  domain.  However,  no 
discussion  is  given  regarding  the  essential  features  governing  the  correct 
response.  Ider  &  Amirouche  [16]  also  develop  an  algorithm  for  multibody 
systems  using  assumed  modes  and  Kane’s  equations.  Their  fonnulation  is 
tailored  for  structures  with  variable  cross-section  beam  elements.  They 
also  consider  the  spin-up  problem  in  their  numerical  examples. 

Taking  an  analytical  approach,  Silverberg  &  Park  [17]  explore 
contributions  of  Coriolis  and  centrifugal  forcing  terms  in  the  response  of 
maneuvering  spacecraft.  Through  development  of  stiffness  operators,  they 
compare  natural  frequencies  of  a  spinning  beam  achieved  by  linearization 
about  both  static  and  dynamic  equilibriums.  They  show  that  linearization 
about  the  dynamic  equilibrium  (same  as  including  geometric  nonlinearities) 
has  an  important  effect  when  certain  nondimensional  spin  and  material 
parameters  are  exceeded. 
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Chapter  2 

Flexible  Body  Formulation 


This  chapter  develops  the  equations  of  motion  for  a  rigid  body  with 
attached  flexible  appendage,  without  articulation,  using  the  principle  of 
virtual  work.  The  governing  equations  are  consistently  derived,  so  that  all 
terms  are  retained. 

Since  the  virtual  work  principle  is  the  basis  for  the  finite  element 
formulation,  it  is  natural  and  consistent  to  discretize  the  flexible  appendage 
using  the  finite  element  method.  In  this  chapter,  lumped  mass 
(3  DOFs/node)  and  lumped  mass/inertia  (6  DOFs/node)  assumptions  are 
employed  in  the  treatment  of  the  mass  distribution  of  the  flexible 
appendage.  Lumped  masses  are  located  at  the  nodes  resulting  from  finite 
element  discretization  of  the  appendage  stiffness.  These  assumptions  lead 
to  equations  of  motion  which  can  be  expressed  in  a  convenient  matrix 
form.  The  exact  governing  equations  could  also  be  fully  developed  using 
finite  element  method,  as  has  been  done  in  the  literature  [14,  15,  18]. 

The  assumption  of  rigid  central  body  allows  governing  equations  of 
motion  to  be  derived  with  respect  to  the  body  fixed  frame.  Alternatively, 
consideration  of  a  free  floating  deformable  body  requires  the  application  of 
conservation  of  linear  and  angular  momentum  [18],  or  other  corotating 
frames  [19,  20]. 

Equations  derived  in  this  chapter  are  quite  general  in  nature,  and  can 
be  applied  to  a  wide  class  of  problems.  Many  satellites,  as  well  as  the  space 
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shuttle,  can  be  approximated  by  the  rigid  central  body  assumption.  The 
only  restriction  on  the  flexible  appendage  is  the  fixed,  or  cantilevered 
boundary  condition,  relative  to  the  rigid  body.  Example  problems  using 
the  lumped  mass/inertia  formulation  are  considered  in  chapter  5,  where  the 
flexible  domain  is  discretized  using  the  beam  finite  elements  derived  in 
chapter  3. 


Figure  2.1.  Rigid  body  with  flexible  appendage. 
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where 


Ti  =  inertial  reference  frame 
f  b  =  body  fixed  reference  frame 
Vr  =  rigid  body  domain 

Vp  =  flexible  body  domain,  reference  configuration 
^R(t)  =  material  coordinate  of  rigid  particle  with  respect  to 
inertial  frame 

^F(t)  =  material  coordinate  of  flexible  particle  with  respect  to 
inertial  frame 

R(t)  =  inertial  position  of  body  frame 
•*R  =  rigid  particle  position  with  respect  to  body  frame 
rF(t)  =  reference  flexible  particle  position  with  respect  to  body 
frame 

c  =  rigid  c.m.  with  respect  to  body  frame 
s(t)  =  vehicle  c.m.  with  respect  to  body  frame 
T|(t)  =  displacement  of  flexible  particle  due  to  deformation 


Figure  2.1  shows  a  rigid  body  with  attached  flexible  appendage  in 
inertial  space.  The  form  of  the  appendage  is  arbitrary,  as  suggested  by  the 
figure.  In  practice,  the  body  fixed  reference  frame  is  located  as  a  matter 
of  convenience,  and  is  not  normally  coincident  with  the  rigid  body  c.m.  or 
the  vehicle  c.m. 


2.2  Principle  of  Virtual  Work 

The  principle  of  virtual  work  is  a  statement  that  for  a  body  in 
equilibrium  under  the  action  of  prescribed  body  and  surface  forces  the 
work  done  by  these  forces  through  a  kinematically  admissible  displacement 
is  equal  to  the  change  in  internal  virtual  work.  In  combination  with 
D'Alembert's  principle,  the  virtual  work  principle  can  be  expressed  as 
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SWext  =  (  -  P^)  dV  +  t  dS  =  (  5e:a  dV  =  SWim 

Js, 


where,  in  addition  to  the  quantities  previously  defined, 

f  =  body  force  (force/volume) 
t  =  surface  traction  applied  over  Sa 

p  =  density 

e  =  strain 

a  =  stress 


(2.2.1) 


The  virtual  work  statement  is  evaluated  as  the  sum  of  two  parts; 
integration  over  the  rigid  central  body  and  integration  over  the  flexible 
appendage. 


2.3  YgctriX-NQtgtiQ.n 

In  forming  the  virtual  work  expression,  several  vectors  must  be 
defined.  When  the  dynamic  system  involves  more  than  one  reference 
frame,  as  is  the  case  for  the  vehicle  of  figure  2.1,  it  is  helpful  to  use  a 
notation  which  explicitly  identifies  the  frame  in  which  vector  components 
are  expressed.  Toward  this  end,  the  vectrix  notation  [21]  is  used.  A  vector 
can  be  written  as  the  multiplication  of  two  column  matrices:  one 
containing  the  vector  components,  the  other  the  frame  unit  directions.  For 
example,  an  arbitrary  vector  v  can  be  expressed  in  some  reference  frame 
f  a»  whose  basis  vectors  are  ai,  32,  is,  as 

v  =  ViSi  +  V2a2  +  V3a3  =  FaY=vT’Fa  (2.3.1) 
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where 


SJ  =  [  Vi  V2  V3  ] 

f* _  r  ^  ^ 

a  -  [  ai  32  33  J 

Differentiation  of  vectors  involving  multiple  reference  frames 
introduces  the  vector  cross  product.  It  is  convenient  to  express  this 
operation  as  a  matrix  multiplication,  in  conjunction  with  vetrices.  The 
cross  product  of  arbitrary  vectors  u  and  v,  both  expressed  in  Tn,  is  given 
by 

T  T 

U  X  V  =  u'TfaXfa  V  =  ra  (u'' v)  (2.3.2) 

where  v  is  the  same  as  above  and  the  components  of  u  have  been  formed 
into  the  skew  symmetric  matrix  given  by 

0  -U3  U2 

u^  =  U3  0  -Ui 

.  -U2  ui  0  J  (2.3.3) 

The  inertial  time  derivative  of  the  frame  is  also  an  important 
relationship  and  is  given  by 

f b  =  b  X  (2.3.4) 
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From  figure  2.1,  the  inertial  positions  of  rigid  and  flexible  material 
particles  are  given  by  the  vectors 


=  R  +  Tr 
=  R  +  rp  +  T1 


(2.4.1) 


(2.4.2) 


Inertial  accelerations  are  obtained  through  time  differentiation  of  the 

position  vectors.  Let  — { )  denote  time  differentiation  in  the  inertial  frame, 

dt 

and  ( )  indicate  time  differentiation  in  the  body  fixed  frame. 
Differentiation  of  equation  (2.4.1)  gives 


Pr  =  d.[dR 
^  dt  L  dt  dt  - 


(2.4.3) 


where  it  is  convenient  to  let  u  =  the  inertial  velocity  of  the  body  fixed 

frame.  Now  expressing  components  of  u  and  in  the  body  fixed  frame 

dt 

allows  the  derivative  to  be  written  as 


dt  L  dt  J 


(2.4.4) 


Application  of  the  chain  mle  and  using  the  vectrix  notation  gives  the  final 
form  of  the  inertial  acceleration  of  an  arbitrary  rigid  material  particle. 
The  acceleration  of  an  arbitrary  flexible  particle  follows  in  an  analogous 


manner. 


•  *  T  r  1  T 

=  Fb  [u  +  (0^  U  +  IR  +  CO^  01^  ngj  =  Sr 


(2.4.5) 
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4f  =  U  +  co^  u  +  (b^  (rp  +  n)  +  fj  +  (0^  {r£  +  T]J  +  5]  =  f  §E 

(2.4.6) 


Virtual  displacements  are  obtained  through  variation  of  the  position 
vectors,  and  can  be  expressed  as 

S^R  =  f [Sx  +  bQ""  rj  (2.4.7) 

=  Tb  [bx  +  bri  +  ba^  (rp  +  gj]  (2.4.8) 


where 


IR,  £F,  a  = 


U  = 

ca  = 

bx  = 
bQ  = 

bn  = 
= 


components  of  tr,  rp,  t|,  expressed  in  the  body  fixed 
frame 

components  of  the  inertial  velocity  of  the  origin  of  the 
body  fixed  frame 

components  of  the  inertial  angular  velocity  of  the  body 
fixed  frame 

components  of  the  variation  of  inertial  position  vector  R 
components  of  angular  variation  which  arises  as  a 
consequence  of  rotation  of  the  body  frame  with  respect  to 
the  inertial  frame 

components  of  the  variation  of  the  relative  displacement 
vector  r| 

b  1  b2  b3  J 


All  components  defined  above  are  expressed  in  the  body  fixed  frame, 
Fb-  Note  that  the  acceleration  of  a  flexible  particle  has  a  high  degree  of 
coupling  between  rigid  body  motions  and  flexible  deformations. 
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The  governing  equations  of  motion  are  obtained  by  substituting 
equations  (2.4.5)-(2.4.8)  into  equation  (2,2.1).  Evaluation  of  equation 
(2.2.1)  produces  terms  which  can  be  grouped  according  to  the  variation 
(5x,  5Q.,  or  bn)  multiplying  each.  Because  of  the  arbitrariness  of  the 
variations,  each  group  must  independently  be  equal  to  zero.  Thus  the 
following  three  sets  of  governing  equations  are  obtained. 

Recall  that  no  restrictions  have  been  imposed  regarding  the  physical 
shape  of  the  flexible  appendage,  or  on  the  form  of  the  relative 
displacements  n.. 
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(0.^  IfAI  =  I  p(lF  +  Tlf  (0^  CO^  (rp  +  T]j  dV 

Ap  (2.5.12) 

along  with  the  following  definitions: 

F  =  components  of  total  force  acting  through  body  frame 
origin,  due  to  body  forces  (JEb)  and  surface 
tractions  (Ft),  expressed  in  the  body  fixed  frame 
T  =  components  of  total  torque  acting  about  body  frame 
origin,  due  to  body  forces  (Tb)  and  surface 
tractions  (It),  expressed  in  the  body  fixed  frame 
I  =  instantaneous  inertia  matrix  of  the  vehicle  about  the  body 
fixed  frame  due  to  rigid  body  (Ir)  and  flexible  body  (Ip), 
expressed  in  body  fixed  frame 
fp  =  components  of  body  forces  acting  on  flexible  body, 
expressed  in  body  frame 

tp  =  components  of  surface  tractions  acting  on  flexible  body, 
expressed  in  body  frame 
m  =  total  vehicle  mass  (mR  +  mp) 

Equations  (2.5.1)-(2.5.3)  are  an  exact  set  of  equations  governing  the 
idealized  vehicle  of  figure  2.1.  Discretization  of  the  flexible  domain  stems 
from  these  equations.  Two  lumped  parameter  approximations  will  be 
considered:  (a)  discretization  of  the  flexible  domain  into  a  collection  of 
point  masses  (no  rotatory  inertia)  interconnected  by  massless  springs,  and 
(b)  extension  of  the  previous  model  to  include  rotatory  inertia. 

In  practice,  mass  is  concentrated  at  locations  corresponding  to  finite 
element  nodes,  which  allows  the  use  of  the  finite  element  stiffness  matrix. 
For  3-D  finite  elements,  condensation  technique  must  be  used  to  make  the 
stiffness  matrix  compatible  with  approximation  (a). 
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Finite  element  discretization  can  also  be  consistently  applied  to  all 
integrals  in  Equations  (2.5.1)-(2.5.3).  Derivation  is  straightforward,  but 
the  solution  procedure  is  more  involved  than  lumped  mass  assumptions, 
and  is  therefore  not  considered  here.  This  approach  has  been  used  to 
derive  equations  of  motion  and  demonstrated  through  simple  problems 
using  quadratic  beam  elements  [18], 


2.6  Lumped  Mass  Assumption 


ith  mass  particle 


2Ti  =  gti- 

Figure  2.2.  Lumped  mass  discretization  of  the  flexible  appendage 


Ui 

Vi 

LWij 


Schematic  representation  of  this  assumption  is  shown  in  figure  2.2, 
where  the  solid  line  represents  the  surface  of  the  flexible  appendage.  The 
flexible  domain  is  modelled  with  a  finite  number  of  point  mass  particles, 
connected  by  massless  springs.  Mathematically,  this  is  a  straightforward 
process  whereby  the  integrals  in  equations  (2.5.1)-(2.5.3)  over  the  flexible 
domain  are  replaced  by  summations  over  the  number  of  mass  particles. 
The  location  of  lumped  masses  correspond  to  the  nodal  locations  of  the 
fmite  element  discretization  of  the  flexible  domain. 
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N 


/to’  i  n)p  dv-^x/toi’i’ 

i=l 


where 


(2.6.1) 


Qi  =  ith  nodal  displacements  (translations  only),  expressed  in 
the  body  fixed  frame 
mi  =  corresponding  lumped  mass 

The  motion  of  a  lumped  mass  is  fully  described  by  three  translations. 
For  3-D  finite  element  models  of  the  flexible  domain  which  include 
rotations  as  nodal  DOFs,  condensation  technique  must  be  used  to  remove 
rotational  DOFs  from  the  stiffness  matrix.  Consistent  with  the  lumped 
mass  approximation,  external  loads  on  the  flexible  body  are  also  Tumped’ 
at  the  nodes.  The  exact  governing  equations  (2.5.1)-(2.5.3)  are  simplified 
through  the  lumped  mass  assumption  to  give 

N  N 

E  =  mil  +  mm.^  U  +  ma^  £  +  mm.^  £  +  2'^  mi  ^  ^  mi  ^ 

i=l  i=l  (2.6.2) 

N 

T  =  m£^  u  +  m£^  (O^U  +  I(a+oi^Im.  +  2^  mi  (r^  +  qjj^  ^ 

i=l 


N 

+  5^  mi(ri  +  q^^ 
i=l 


(2.6.3) 


fi  =  mill  +  mim.^  u  +  mifll^  (li  +  Si)  +  miia.^  (ri  +  ^  +  2mim.^  ^ 


N 

-I-  mi^  +  X  9i  ,2,...,N 

j=i  (2.6.4) 
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where 


li  =  reference  position  of  ith  node,  components  in  the  body 
fixed  frame 

£  =  lumped  nodal  forces 

Kij  =  assembled  stiffness  matrix  (rotation  condensed  out) 

Equations  (2.6.2)-(2.6.4)  are  valid  for  nonlinear  flexible  systems  by 
interpretation  of  the  stiffness  matrix  as  the  tangent  stiffness  matrix.  An 
alternative  derivation  of  the  lumped  mass  equations  of  motion  are  provided 
in  [22]. 


2.6.1  Equations  of  Motion 

Equations  (2.6.2)-(2.6.4)  are  the  lumped  mass  equations  of  motion, 
comprising  (6+3N)  scalar  equations.  These  equations  can  be  recast  into  a 
single  matrix  equation  of  motion  which  can  be  numerically  integrated. 

Mil  +  KU  =  R  (2.6.5) 


It  is  natural  to  partition  the  matrix  equation  into  rigid  and  flexible  body 
contributions  and  rewrite  equation  (2.6.5)  as 
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(2.6.7) 
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(2.6.9) 

(2.6.10) 

(2.6.11) 

(2.6.12) 
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Rrf  = 
(6x1) 


N 

-2^  mico''  ^ 

i=l 

N 

-2]^  mi(ri  +  qif  (O'"  ^ 
i=l 


Rf 

(3Nxl) 


-  mico^  u  -  miO)^  (fj^  +  qj  -  2miOi^  ^ 


-  niNfii^  U  -  niNOi^  (lii  +  qw)  -  2mNffi.^  ^ 


(2.6.13) 


(2.6.14) 


Note  that  the  mass  matrix  M  is  symmetric,  positive  definite;  the 
stiffness  matrix  K  is  symmetric,  positive  semi-definite  and  allows  rigid 
body  motion.  Numerical  solution  can  be  obtained  in  a  number  of  ways. 
For  linear  systems,  transformation  can  be  made  to  modal  space,  which 
allows  truncation  of  both  system  size  and  high  frequency  modes.  Modal 
reduction  is  generally  not  possible  for  nonlinear  flexible  systems.  Thus 
direct  integration  is  preferred  in  the  present  context.  Also  the  effect  of 
various  terms  can  be  more  readily  assessed  in  physical  space. 


The  methods  available  for  direct  integration  of  equation  (2.6.5)  can 
be  classified  into  explicit  and  implicit  schemes.  Explicit  schemes  are 
conditionally  stable  and  require  very  small  time  steps  to  integrate  the 
highest  frequencies  accurately  [23,  24].  Implicit  schemes  are  advantageous 
because  they  are  unconditionally  stable  and  the  step  size  can  be  chosen  on 
the  basis  of  accuracy  only.  This  generally  allows  a  much  larger  step  size 
than  would  be  required  by  explicit  schemes.  Because  of  the  relaxed 
integration  step  size  afforded  by  unconditionally  stable  implicit  schemes, 
the  Newmark  integration  method  is  implemented  in  the  examples  of 
chapter  5.  Details  of  the  Newmark  scheme  are  outlined  in  Appendix  B. 
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2.7  Lumped  Mass/Inertia  Assumption 


Figure  2.3.  Lumped  mass/inertia  discretization  of  flexible  appendage. 


As  in  the  previous  derivation,  integrals  over  the  flexible  domain  are 
replaced  by  summations  over  a  finite  number  of  nodes.  The  nodes  are  now 
treated  as  small  rigid  bodies;  they  are  no  longer  mass  particles,  but  have 
both  mass  and  inertia.  Six  quantities  are  required  to  describe  the  motion  of 
each  node. 

Massless  springs  connect  the  nodal  bodies  just  as  in  section  2.6.  In 
practice,  the  finite  element  stiffness  matrix  provides  information  regarding 
interconnecting  forces.  The  lumped  mass/inertia  formulation,  as 
constrasted  with  the  lumped  mass  formulation,  has  the  advantage  that  finite 
element  DOFs  are  used  directly;  no  condensation  is  required  in  the 
numerical  solution. 


37 


2.7.1  Kinematic  Description 

undeflected  deflected 

configuration  configuration 


Figure  2.4.  Relative  displacement  of  ith  nodal  body. 


Consider  the  rigid  nodal  body  shown  in  figure  2.4  with  an  arbitrary 
material  particle  labelled  “a”  in  the  undeflected  configuration,  “A”  in  the 
deflected  configuration.  As  measured  with  respect  to  the  body  fixed 
frame,  the  nodal  body  undergoes  infinitesmal  translation  and  rotation  in 
moving  to  the  deflected  configuration.  The  location  of  an  arbitrary 
material  particle  in  the  undeflected  configuration  is  given  by 

rF  =  rFo  +  >-  (2.7.1) 

where  rpo  is  the  location  of  the  nodal  body  reference  frame  in  the 
undeflected  configuration.  The  nodal  body  reference  frame  is  aligned  with 
the  body  fixed  reference  frame  in  the  undeflected  configuration. 
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In  the  deflected  configuration,  the  location  of  A  is  given  by 

Rf  =  Rfo  +  A  =  rpo  +  qi  + A  (2.7.2) 

where  Rfo  is  the  location  of  the  nodal  body  reference  frame  in  the 
deflected  configuration.  The  relative  displacement  undergone  by  a 
material  particle  in  moving  from  a  to  A  is  then 

Ti  =  RF-rF  =  qt  +  A-X  {1.13) 


Rotation  of  the  nodal  body  reference  frame  can  be  expressed  by  the 
skew  symmetric  infinitesmal  rotation  matrix  so  that 


where 


-9z  0y 

1  -Ox 
ex  1 


(2.7.4) 


Substitution  of  equation  (2.7.4)  into  (2.7.3)  leads  to  expression  of  the 
relative  displacement  as 

T|  =  qt  +  -  X  =  qt  +  (C  -  I)X  (2.7.5) 


where  I  is  the  identity  matrix.  Now  since  all  vectors  are  expressed  in  the 
body  fixed  frame,  the  component  notation  is  adopted.  The  matrix  (C  - 1) 
above  can  be  compared  with  the  skew  symmetric  matrix  (associated  with 
vector  cross  products)  introduced  in  section  2.3.  Define 
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so  that  =  (C  - 1).  Finally,  the  relative  displacement  of  a  material 
particle,  assuming  infinitesmal  rotation,  can  be  written  as 

Ti  =  qt  +  qr''2.  (2.7.6) 


2.7.2  Equations  of  Motion 

The  equations  of  motion  follow  from  substitution  of  equation  (2.7.6) 
into  equations  (2.5.1)-(2.5.3)  and  evaluating.  Reduction  of  the  exact 
equations  is  complicated  somewhat  by  the  form  of  the  relative  displacement 
given  by  equation  (2.7.6).  Some  integrals  produce  higher-order  terms 
which  cannot  be  given  simple  physical  interpretation.  These  higher-order 
terms  are  ignored.  Details  of  integrations  are  not  presented,  and  the 
lumped  mass/inertia  equations  of  motion  in  matrix  form  are  given  by 


'Mrr 
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’  Ur  ' 
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■  Ur  ■ 

Rr  +  Rrf 
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Mff  . 

.  tip  . 

0 

1 _ 

.  Up 

Rf 

where  all  partitions  are  explicitly  defined  below.  Note  that  the  rigid-rigid 
partition  is  the  same  as  in  section  2.6.  Kff  is  the  full  stiffness  matrix 
produced  by  the  finite  element  method.  The  tangent  stiffiiess  matrix  is 
used  for  the  solution  of  problems  with  nonlinear  flexibility. 
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(2.7.8) 
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E  -  mco^  u  -  £ 

Kr  = 

(6x1)  L  I  -  nic^  co^  u  -  (0^  I  w 


Rrf  - 
(6x1) 
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i=l 


N  N 

-2^  mi{ri  +  ^  +  ^  li^ 

i=l  i=l 


(2.7.13) 


(2.7.14) 


Rf 

(6Nxl) 


fi  -  mico^  u  -  mi(0^  (n  +  £yj  -  2mi(0^  ^  +  h.o.t. 

ti  +  h.o.t. 


-  niNiO''  u  -  niNm.^  cfl.^  (fN  +  g[tN)  -  2mN0i^  ^  +  h.o.t. 

tN  +  h.o.t. 


(2.7.15) 


2.8  Extending  Symbolic  Rigid  Body  Codes 

Implementation  of  the  flexible  body  formulation  can  take  advantage 
of  available  symbolic  rigid  body  software.  These  software  packages 
produce  FORTRAN  coding  of  the  equations  of  motion  of  a  user  specified 
rigid  multibody  system.  Some  examples  of  symbolic  manipulation  rigid 
multibody  software  include  SD/FAST  [25],  AUTOSIM  [26],  and 
AUTOLEV  [27].  The  use  of  symbolic  rigid  body  codes  allows  the  analyst 
to  concentrate  on  a  smaller  set  of  ‘hand  derived’  equations  addressing  the 
flexible  domain  [28,  29]. 

To  show  how  the  rigid  body  subroutines,  generated  by  any  one  of 
the  above  programs,  can  be  used  in  the  flexible  body  implementation,  the 
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partitioned  equations  of  motion  (2.6.6)  or  (2.7.7)  are  rewritten  as  a  pair  of 
matrix  equations 


Mfr  Ur  +  MpF  tip  +  KpF  Up  =  Rp  (2.8.1b) 

where  it  is  noted  that  Mrr  and  Rr  refer  to  the  current  configuration  of  the 
vehicle,  i.e  the  rigid  body  +  flexible  appendage  is  assumed  rigid  in  the 
current  configuration. 

The  rigid  body  code  produces  a  set  of  subroutines  to  solve  the  rigid 
equations  of  motion 

Mrr  Ur  =  Rr  (2.8.2) 


Mrr  Ur  +  Mrf  Up  =  Rr 


+  R 


RF 


(2.8.1a) 


2,8.1  For  Second  Order  Integration  Schemes 

When  a  second  order  integration  scheme  is  used  in  the  solution  of  the 
vehicle  equations  of  motion,  the  rigid  subroutine  is  necessary  only  in  the 
calculation  of  Mrr  and  Rr.  The  current  configuration  vehicle  c.m.  and 
inertia  matrix  are  provided  as  inputs  to  the  subroutine.  If  the  vehicle 
undergoes  large  flexible  displacements,  the  configuration  must  be  updated 
at  each  integration  step. 
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2.8.2  For  First  Order  Integration  Schemes 

A  sketch  of  the  implementation  of  a  Runge-Kutta  integration  scheme 
in  the  solution  of  the  equations  of  motion  is  offered  in  Appendix  B.  It  is 
shown  that  the  computational  effort  is  concentrated  on  the  evaluation  of 
derivatives,  in  order  to  employ  the  formula  given  by  equation  (B.4.3).  In 
section  B.5,  equations  (2.8.1)  are  rearranged  to  yield 

[Mrr  -  Mrf  Mff’^  Mfr]  Ur  =  Rr  +  Rrf  +  •  •  •  (B.5.4) 


and 


Uf  =  Mff'^ 


(B.5.3) 


To  obtain  the  derivatives  Ur  and  Uf,  these  two  equations  must  be  solved  in 
the  order  given.  Examination  of  equation  (B.5.4)  indicates  that  the  code 
produced  to  solve  equation  (2.8.2)  can  be  used  in  the  flexible  context  by 
modification  to  Mrr  and  Rr.  Simulations  using  first  order  integration 
schemes  can  take  advantage  of  the  code  produced  by  symbolic  software  to  a 
larger  extent  than  simulations  using  second  order  schemes. 
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This  chapter  presents  the  derivation  of  two  nonlinear  beam  finite 
elements.  Isotropic,  prismatic  beam  elements  are  allowed  to  stretch,  bend 
in  two  planes,  and  undergo  St.  Venant  torsion.  For  the  Timoshenko  beam, 
shearing  in  two  planes  is  also  allowed.  To  derive  the  beam  finite  elements, 
assume  that  the  flexible  domain  of  chapter  2  (see  figure  2.1  and  equations 
(2.5.1)-(2.5.3))  is  a  beam,  which  allows  explicit  statement  of  kinematic 
assumptions.  Bemoulli-Euler  kinematic  assumptions  comprise  the 
“engineering  theory  of  beams.”  A  distinction  is  made  between  Bemoulli- 
Euler  and  Rayleigh  beam  theories  in  dynamics  (Bemoulli-Euler  ignores 
rotatory  inertia).  The  Timoshenko  kinematic  assumptions  lead  to  a  beam 
theory  which  includes  the  effects  of  rotatory  inertia  and  shear  strain  within 
the  beam. 

The  development  proceeds  from  the  3-D  statement  of  the  principle 
of  virtual  work.  Kinematic  assumptions  are  explicitly  introduced,  and  the 
work  expression  is  integrated  across  the  beam  cross-section  area  for 
reduction  to  a  1-D  theory. 

The  beam  elements  derived  are  generally  known  as  ‘isoparametric’ 
elements;  isoparametric  meaning  ‘same  parameters’.  For  the  current 
discussion,  let  element  geometry  be  interpolated  from  nodal  values  by 
using  shape  functions.  Let  element  displacements  be  interpolated  from 
nodal  DOFs  by  using  interpolation  functions.  Strictly  speaking. 


isoparametric  formulations  employ  the  same  function  for  shape  and 
interpolation.  It  should  be  noted  that  although  the  Cl  element  to  follow 
does  not  adhere  to  this  rigorous  definition  (it  uses  both  linear  and  cubic 
interpolation  functions),  it  may  be  loosely  referred  to  as  isoparametric. 
The  CO  element  is  a  true  isoparametric  element. 


3.1  Principle  of  Virtual  Work 


The  principle  of  virtual  work  has  been  previously  stated  in  section 
2.2,  where  it  was  applied  to  a  vehicle  composed  of  rigid  central  body  and 
attached  flexible  appendage.  Three  coupled  sets  of  integral  equations  were 
derived  which  govern  the  vehicle  motion.  Equations  (2.5.1)  and  (2.5.2) 
govern  the  translation  and  rotation  of  the  the  rigid  body.  Equation  (2.5.3) 
governs  the  deflections  of  the  flexible  appendage,  relative  to  the  body  fixed 
frame.  For  the  purpose  of  this  derivation,  consider  the  rigid  body  to  be 
fixed  in  inertial  space,  i.e.  fli  =  il  =  0.  Thus  surviving  terms  give  the 
virtual  work  expression  which  governs  the  deflections  of  the  beam  relative 
to  a  body  fixed  reference  frame. 


+ 


5]!^^  dS  = 


Sq^^iip  dV  + 


5e:adV 


where 


(3.1.1) 


H  =  components  of  flexible  displacement,  with  respect  to  the 
body  fixed  frame 

fp  =  force/unit  volume,  with  respect  to  the  body  fixed  frame 
iF  =  surface  traction  applied  over  Sof*  with  respect  to  the  body 
fixed  frame 
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p  =  density 

Vp  =  reference  configuration  of  the  flexible  body 

e  =  strain 

o  =  stress 

3.2  Beam  Theory  Preliminaries 

X,  u 


Figure  3.1.  3-D  Beam  Coordinate  System. 

Consider  the  beam  shown  in  figure  3.1,  with  coordinate  axes  as 
shown.  It  is  assumed  that  stresses  Gy  and  Gz  are  small  compared  to  Gx.  St. 
Venant  torsion  theory  is  incorporated  into  the  displacement  field  (cross- 
sections  undergo  zero  warping,  and  there  is  no  distortion  of  the  cross- 
section).  The  displacement  field  associated  with  the  Bemoulli-Euler 
kinematics  are 

u(x,y,z)  =  Uo(x)  -  yvo,x(x)  -  zwo,x{x} 

v(x,z)  =  Vo(x)  -  z9x(x)  (3.2.1) 

w{x,y)  =  Wo(x)  +  y0x{x) 

and  can  be  interpreted  geometrically  to  mean  that  cross-sections  remain 
perpendicular  to  the  neutral  axis  during  deformation.  A  representative 
construction  is  shown  in  figure  3.2. 
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Figure  3.2.  Kinematics  of  Bemoulli-Euler  Beam  Theory. 


In  the  Timoshenko  beam  theory,  cross-sections  initially 
perpendicular  to  the  neutral  axis  of  the  beam  remain  plane  but  not 
perpendicular  to  the  neutral  axis  during  deformation.  This  kinematic 
assumption  is  shown  geometrically  in  figure  3.3. 


Figure  3.3.  Kinematics  of  Timoshenko  Beam  Theory. 
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The  displacement  field  associated  with  the  kinematics  of  Timoshenko 
beam  theory  are 

u(x,y,z)  =  Uo(x)  -  y0z(x)  +  z0y(x) 

v{x,z)  =  Vo(x)  -  z0x(x)  (3  2.2) 

w(x,y)  =  Wo(x)  +  y0x(x) 

where  zero  subscripts  indicate  displacement  of  the  neutral  axis  and  0’s  are 
rotations  about  the  subscripted  axes. 

Note  the  difference  in  sign  of  the  components  in  the  displacement 
field  u(x,y,z)  of  Bemoulli-Euler  and  Timoshenko  kinematics.  Positive 
rotation  associated  with  w,x  disagrees  with  the  right  hand  rule,  and  as  a 
consequence,  also  disagrees  with  the  right  handed  coordinate  system  used  in 
the  derivation  of  the  vehicle  equations  in  chapter  2.  Therefore,  in  order  to 
use  the  Bemoulli-Euler  (Cl)  element  in  the  dynamic  simulation,  a 
transformation  must  be  made  so  that  the  rotation  is  consistent  with  a  right 
hand  coordinate  system. 


The  restriction  of  zero  deformation  of  the  cross-section  during 
torsion  implies  that  Yyz  =  0.  The  assumptions  associated  with  both 
kinematic  models  allow  reduction  of  the  3-D  linear  elastic,  isotropic  stress- 
strain  relations  to  the  simple  result 


(3.2.3) 
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3.3  C*  Formulation 

The  formulation  encompasses  both  the  Bemoulli-Euler  and 
Rayleigh  beam  theories.  The  subtle  distinction  is  the  inclusion  of  rotatory 
inertia  in  Rayleigh’s  equations  governing  the  motion  of  a  beam  [30].  It  will 
be  shown  that  consistent  derivation  using  the  virtual  work  principle 
produces  the  Rayleigh  theory.  By  convention,  Bemoulli-Euler  theory 
excludes  rotatory  inertia. 

From  small  deflection  theory,  the  non-zero  strains  associated  with 
the  displacement  field  of  equation  (3.2.1)  are 

ex  =  U.X  =  Uo.x(x)  -  yvo,xx{x)  -  zwo^x(x) 

7xy  =  U  y  +  V,x  =  -  zex,x(x)  (3  3  1) 

7zx  =  W.X  +  u^  =  yex,x(x) 

Substituting  equations  (3.3.1)  and  (3.2.1)  and  constitutive  relations 
(3.2.3)  into  the  right  hand  side  of  (3.1.1)  gives  the  internal  virtual  work 
expression  in  the  volume  integral  form 

j  |p[(5uo  -  y5vo,x  -  z5wo J  (iio  -  yvo,x  -  zwoj  +  (5vo  -  zbBx)  (vq  -  zSx) 

+  (5wo  +  ySGx)  (wo  +  y^x)] 

+  E(5uo,x  -  y5Vo,xx  -  z5Wo,xx)(Uo,x  -  yVo,xx  '  ZWq^ix) 

+  G(-z50xJ  (-zGx^)  -h  G{y50x,x)(yex,x))  dV  =  RHS  (3.3.2) 

where  the  domain  is  understood  to  be  the  flexible  volume.  The  expression 
(3.3.2)  can  be  integrated  through  the  beam  cross-section  to  give 
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(m(5UoUo  +  5VoVo  +  5WoWo)  +  (imy  +  Imz)50x0x 

Imy5Wo,xWo,x  ImzS^o,x^o,x  EASUo,xUo,x 
+  EIy5Wo,xx Wo,xx  +  EIz6Vo,xxVo,xx  +  ^(ly  Iz)S0x,x0x,x)dx  =  RHS  (3.3.3) 

where  cross-section  principal  axes  are  assumed  to  be  aligned  with  element 
coordinate  axes.  Constants  appearing  in  equation  (3.3.3)  are  defined  as 
follows: 


m  =  massAength 
ly,  Iz  =  area  moment  of  inertia 
Imy»  Imz  =  mass  moment  of  inertia 

The  integral  equation  (3.3.3)  appearing  above  is  to  be  evaluated  over 
the  length  of  the  flexible  domain.  The  flexible  domain  is  decomposed  into 
finite  elements  and  equation  (3.3.3)  becomes  a  summation  of  integrals  over 
the  individual  element  domains.  The  continuity  of  assumed  displacement 
functions  is  dictated  by  the  terms  appearing  in  the  integrand.  Since  second 
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derivatives  of  displacement  exist  in  the  element  stiffness  integral,  inter¬ 
element  continuity  requires  that  the  displacement  fields  w-,  and  v  be 
continuous.  The  axial  displacement  and  twist  require  only  continuity. 
(The  superscript  indicates  the  derivative  through  w'hich  the  function  or 
field  is  continuous.)  The  continuity  condition  is  the  integrability  of 
equation  (3.3.3)  over  the  entire  domain. 

3.3.1  Discretization 

^=-1  i=+i 

Node  Node 

1  2 

Figure  3.4.  Local  coordinates  for  2  node  beam  element. 


The  two  node  element  is  defined  by  the  local  coordinate  system 
shown  in  figure  3.4.  The  element  has  six  DOFs  at  each  node.  The 
necessary  interpolation  functions  are  defined  in  terms  of  the  local 
coordinate  The  linear  interpolation  functions  are  given  by 


N?  =  I(l-5) 

N5  =  i(i  +  y 


(3.3.4) 


For  the  Cl  continuous  displacement  fields  w©  and  Vq,  the  cubic 
Hermitian  interpolation  functions  are  used  and  are  defined  as 

Nl  =  i(l-y(2-5-4^) 

4 
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Ni=Mi-4ni+ii 

4 

N]  =  i(i  +  y(2  +  4-52) 

4 

4 


Element  nodal  DOFs  are  arranged  as 


II 

where 

^  =  [  Uol 

Vol  Woi 

0x1 

W0.XI  Vo,xl  ] 

92'^  =  [  Uo2 

Vo2  Wo2 

0x2 

Wo^2  Vo,x2  ] 

Also  define 

=  [  Uo  Vo  Wo  0x  Wo.x  Vo.x  ] 

~  [  Uo,x  0x,x  'Vo,xx  Vo,xx  ]  ~  [  £ox  ] 

where  the  strains,  are  defined  as  the  work  conjugates  to  the  stress 
resultants  (see  Appendix  A).  These  definitions  allow  a  succinct  way  of 
defining  the  matrices  N.  and  B.,  which  relate  the  element  nodal 
displacements  to  the  vectors  ^  and  £. 

With  the  interpolation  functions  and  nodal  DOFs  explicitly  stated, 
some  explanation  is  necessary  to  avoid  a  pitfall  which  can  be  harmlessly 
overlooked  in  elements.  Nodal  DOFs  consist  of  axial  displacement  and 
transverse  displacements  and  their  derivatives  with  respect  to  x,  while  the 
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interpolation  functions  are  defined  in  local  coordinates.  Care  must  be 
taken  to  insure  correct  application  of  the  Jacobian  transformation. 

To  consider  the  problem  of  interpolation  in  a  consistent  manner, 
note  that  interpolation  is  performed  using  the  nodal  quantities  Wo,^i  and 
For  example. 

Wo  =  nIwoi  +  N2Wo,4i  +  N3W02  +  n1wo.42 

W0.4  =  Ni,^Wol  +  Ni^Wo,4l  +  N3,^Wo2  +  N4.4Wo,^2  (3.3.7) 

Wo,^^  =  Ni,^^Woi  +  N2,44Wo,4i  +  N3,^4Wo2  + 


Note  that  these  expressions  yield  the  expected  results  when  evaluated 
at  the  endpoints,  4  =  ±  1-  In  order  to  express  each  of  these  equations 
properly  in  physical  coordinates,  we  must  consider  how  derivatives 
transform  between  local  and  physical  coordinates.  Since  Ni  =  Ni  (^),  the 
derivative  with  respect  to  x  can  be  written 
dNj  _  dNj  d^ 

dx  d^  dx  -  (33,g) 


where  by  definition,  the  Jacobian  is  J  = 


The  second  derivative  becomes 


d^Ni  _  d^Nj  (d^l^  ^  dNj  d^^ 
dx^  'dxl  d^  dx^ 


Likewise,  displacement  derivatives  can  be  written 

dwp  _  dwp  d^ 
dx  d^  dx 


(3.3.10) 
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(3.3.11) 


where  for  straight  beams,  the  Jacobian  is  constant  and  equal  to  L/2  in  each 
element  (L  is  the  physical  length  of  an  arbitrary  element),  so  that  the 
second  term  of  (3.3.9)  and  (3.3.11)  is  zero.  Substituting  equations  (3.3.8)- 

(3.3.11)  into  equations  (3.3.7),  displacements  and  displacement  derivatives 
can  be  rewritten  in  terms  of  physical  coordinates.  This  gives,  in  matrix 
form. 


Wo  ■ 
Wo,x 

- 1 

‘x 

JNi 

Ni 

Nix 

JNi 

JNix 

Wol 

Wo.xl 

Wo,xx. 

_  nU, 

JNU 

Nlxx 

- 1 

.  s 

z. 

>—> 

Wo2 

.Wo.x2. 

(3.3.12) 


and  demonstrates  correct  differentiation  (with  respect  to  the  physical 
coordinate)  using  local  cubic  interpolation  functions.  A  similar  statement 
can  be  made  for  the  displacement  Vq.  For  displacement  fields  requiring 
only  linear  interpolation  functions,  the  relationship  between 
displacement/displacement  derivatives  and  nodal  DOFs  does  not  involve  the 
Jacobian.  For  instance. 


Uo 

'  N? 

N§  ’ 

Uol 

Uo,X 

.  N?.x 

N?,,. 

.  Uo2  . 

(3.3.13) 


Now  equation  (3.3,3)  can  be  rewritten  in  matrix  form  as 


5d^  Dm  d  +  ^  Dk  £  dx  =  RHS 


(3.3.14) 
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and  finally,  equation  (3.3.14)  can  be  rewritten  in  terms  of  the  nodal  DOFs 
to  give 

6q'^Ma+5q'^Kq=RHS 

where 


M  =  consistent  mass  matrix 
K  =  stiffness  matrix 


and  the  element  level  consistent  mass  matrix  and  material  stiffness  matrix 


are  given  by 
M  =  f  Dm  N  dx 


B^DkBdx 


(3.3.15) 

(3.3.16) 


where 
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and 


EA 

GJ 

Ely 

EIz 

The  matrices  N  and  B.  relate  element  nodal  displacements  to  the 
vectors  d  and  £  and  are  defined  by  the  equations 

d  =  (3.3.17) 

£  =  ^3  (3.3.18) 


These  matrices  are  formed  by  appropriate  placement  of  the  linear  and 
Hermitian  interpolation  functions  and  their  derivatives.  For  the  Cl  element 
formulation,  they  are 

N?  000  0  ON^OOOO  0 

0  N{  0  0  0  JN^  0  0  0  0  JN^ 

0  0  nJ  0  JNJ  0  0  0  0  JN^  0 

N  = 

000  N^O  OOOON^O  0 

0  0  0  0  0  0  0  0 

0  N|  ^  0  0  0  0  ^  0  0  0  JNi  ,  _ 
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and 


N?,x  0  0  0  0 


0  N?,,  0  0  0  0 


0 


B  = 


0  0  0  0  0 

0  0  N\,„  0  JNi,„  0 


0  0  0  0  0 

0  0  n5,„  0  Jn1,„  0 


0  Ni„  0  0  0  JNi,„  0  n5.„  0  0  0  JNj,„ 


The  preceding  equations  are  properly  known  as  the  beam  theory  of 
Rayleigh.  By  convention,  Bemoulli-Euler  theory  sets  the  rotatory  inertias 
Imy  =  Imz  =  0- 

The  above  expressions  for  element  level  consistent  mass  matrix  and 
the  material  stiffness  matrix  contain  local  interpolation  functions 
differentiated  with  respect  to  the  physical  coordinate  x.  Proper 
introduction  of  the  Jacobian  yields  integrals  in  a  form  which  can  be 
numerically  integrated  using  Gaussian  quadrature  (Appendix  F). 
Assembled  matrices  are  obtained  by  summing  over  all  element  matrices. 

3.4  CO  Formulation 

The  CO  formulation  is  derived  based  on  the  kinematic  assumptions  of 
Timoshenko  beam  theory.  Rotation  of  beam  cross-sections  are  independent 
of  transverse  displacements  (recall  figure  3.3).  This  produces  a  beam 
theory  which  has  nonzero  shearing  strains.  This  formulation  is  more 
accurate  when  the  thickness  of  the  beam  is  large  compared  to  its  length, 
and  for  higher  vibration  modes  (wavelength/thickness  — >  small). 
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From  small  deflection  theory,  the  non-zero  strains  associated  with 
the  displacement  field  of  equation  (3.2.2)  are 

Ex  =  U,x  =  Uo,x(x)  -  y0z^{x)  -I-  Z0y.x(x) 

7xy  =  u,y  +  V.X  =  -  0z{x)  Vo.x{x)  -  Z0xji(x)  (3.4.1) 

Yzx  =  W,x  +  U,z  =  Wo.x(x)f  y0x,x(x)  +  0y(x) 

Substituting  equations  (3.4.1)  and  (3.2.2)  and  constitutive  relations 
(3.2.3)  into  the  right  hand  side  of  equation  (3.1.1)  gives  the  internal  virtual 
work  expression  in  the  volume  integral  form 

j  ip[(5uo  -  y60z  +  z50y)  (iio  -  ySz  +  z0y)  +  (6vo  -  z50x)  (vo  -  z6x) 

+  (5wo  +  y50x)  (wo  +  ySx)] 

+  E(5uo,x  -  y50z,x  +  z60y.x){uo.x  -  y0z.x  +  Z0y,x) 

+  G(-50z  +  5vo,x  -  z50x.x)  (-0Z  +  Vo.x  -  Z0x,x) 

+  G(50y  +  5wo^  +  y50x,x)  (0y  +  Wo,x  +  y0x,x)|  dv  =  RHS  (3.4.2) 


where  the  domain  is  understood  to  be  the  flexible  volume.  The  expression 
(3.4.2)  can  be  integrated  through  the  beam  cross-section  to  give 


(m(6UoUo  +  5VoVo  +  SWqWo)  +  (Imy  +  Imz)S0x0x  Imy50y0y  +  ImzS0z0z 

-I-  EA5Uo,xUo,X  +  EIy50y,X0y,X  ■*"  EIz50z,X0Z,X  ^(^y  Iz)S0x,X0X,X 
+  GA8(vo,x  -  0z)(vo,x  -  0z)  +  GA8(wo^  +  0y)(wo,x  +  0y)ldx  =  RHS  (3.4.3) 
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where  cross-section  principal  axes  are  assumed  to  be  aligned  with  element 
coordinate  axes.  Constants  appearing  in  equation  (3.4.3)  are  the  same  as 
defined  in  section  3.3  (see  page  51). 

The  integration  to  be  carried  out  in  equation  (3.4.3)  is  again 
evaluated  as  the  summation  over  the  finite  elements  which  constitute  the 
flexible  domain.  Examination  of  the  integrand  reveals  that  only  C® 
continuity  is  required. 


3.4.1  Discretization 


A  schematic  of  the  finite  element  model  is  shown  in  figure  3.4. 
Because  cross-sectional  rotations  are  independent  of  the  transverse 
displacements,  only  the  linear  interpolation  functions  are  necessary.  In 
local  coordinates  the  interpolation  functions  have  the  form 


N?  =  I(1-|) 

n5  =  1(i  +  4) 


(3.4.4) 


Element  nodal  DOFs  are  arranged  as 

3''  =  [9l^ 


(3.4.5) 


where 

^  =  [  Uol  Vol  Wol  0x1  0yl  0zl  ] 

=  [  Uo2  Vo2  Wo2  0x2  6y2  0z2  ] 
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The  CO  displacements  and  strains  are  given  by 


d  —  [  Uq  Vq  Wq  0x  0y  6z  ] 

~  [  Uo,x  0x,x  0ype  0z,x  ('^o,x  “  0z)  ■*"  ®y)  ] 


[  £xo  “^y  Yxyo  Yxzo  ] 


where  the  additional  terms  in  this  strain  vector,  compared  to  the  previous 
section,  are  the  shear  strains  arising  from  the  Timoshenko  kinematics. 
Forming  these  vectors  from  the  element  nodal  displacements  defines  the 
matrices  N  and  B.  These  matrices  are  composed  of  the  linear  interpolation 
functions  and  their  derivatives. 


(3.4.6) 


£  =  B  q 


(3.4.7) 


Explicitly,  the  N  and  B  matrix  are  given  by 


N?  0  0  0 

0  N?  0  0 

0  0  N?  0 

0  0  0  N? 

0  0  0  0 

0  0  0  0 


0  0  1^5  0 

0  0  0  n5 

0  0  0  0 

0  0  0  0 

N?  0  0  0 

0  N?  0  0 


0  0  0  0 
0  0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  0 
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0  0  0  0  0  n5.x  0  0  0  0  0 

0  0  0  N?,,  0  0  0  0  0  0  0 

0  0  0  0  N^,0  0  0  0  0  N^,0 

a= 

0  0  0  0  0  N?,  0  0  0  0  0  n5.x 

0  N^,x  0  0  0  -N?  0  n5.x  0  0  0  -n5 

_  0  0  N?.x  0  N?  0  0  0  0  0  _ 

Now  equation  (3.4.3)  can  be  rewritten  in  matrix  form 

r 

I  5d'^Dnid  +  5£'^^£dx  =  RHS 

A  (3.4.8) 

and  finally,  equation  (3.4.8)  can  be  rewritten  in  terms  of  the  nodal  DOFs  to 
give 

5q'^M£+6c['^Kq=RHS 

where 


M  =  consistent  mass  matrix 
K  =  stiffness  matrix 


and  the  element  level  consistent  mass  matrix  and  material  stiffness  matrix 


are  given  by 


Dm  N  dx 


B'^DkEdx 


(3.4.9) 

(3.4.10) 
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where 


m 


Dm  — 


m 


m 


(imy  +  Itnz) 


imy 


^mz 


and 


Dj,= 


EA 


GJ 


EL 


EL 


GA 


GA 


The  expression  for  the  material  stiffness  matrix  contains  local 
interpolation  functions  differentiated  with  respect  to  the  physical 
coordinate  x.  Jacobian  transformation  must  be  introduced  in  order  to 
properly  express  the  integrand  in  terms  of  x.  The  integrals  can  then  be 
numerically  integrated  using  Gaussian  quadrature  (Appendix  F). 
Assembled  matrices  are  obtained  by  summing  over  all  element  matrices. 


3.5.  Geometric  Stiffness  Matrix 

For  problems  involving  geometric  nonlinearities,  system  equilibrium 
must  be  satisfied  in  the  deformed  configuration.  In  terms  of  the  virtual 


work  principle,  the  deformed,  or  current  configuration,  is  represented  by 
consideration  of  both  linear  and  nonlinear  strain  terms.  Derivation  of  the 
geometric  stiffness  matrix  is  accomplished  through  incremental 
linearization  of  the  internal  virtual  work  principle  at  the  current 
configuration  of  the  flexible  domain  [31].  All  terms  are  retained  in  the 
consistent  derivation.  The  incremental  virtual  work  statement  is 


A5Wint  =  I  (6 Ae :  Aa  +  5  Ae :  Ot)  dV 

/  Vf 

=  I  6A^(AfF  -  pA^  dV  +  I 

y  vp  Js. 


dV+  I  5A5AtFdS  =  A6Wext 

SoF  (3.5.1) 


where 


Aeij 


2  \  9xj  9xi  / 


Aejj 


pAukj 

2  I  9xi  /  \  axj  I 


cJr  =  total  accumulated  stress 

Afp  =  increment  in  prescribed  body  forces,  components 
expressed  in  body  fixed  frame 

AtF  =  increment  in  prescribed  tractions,  components  expressed 
in  body  fixed  frame 


The  first  term  of  the  incremental  internal  virtual  work  above 
becomes  the  material  stiffness  matrix  previously  obtained.  The  second 
term  leads  to  the  geometric  stiffness  matrix  and  will  be  derived  in  this 
section.  The  nonlinear  strain  increments,  Aeij,  may  be  formed  from  either 
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set  of  kinematic  assumptions  without  great  difficulty.  In  accordance  with 
thin  beam  assumptions  noted  earlier,  Oy  and  Oz  are  small  compared  to  Ox 
and  may  be  neglected.  St.  Venant  torsion  theory  does  not  allow  distortion 
of  the  cross-section,  which  implies  that  Yyz  =  0  and  leads  to  Xyz  =  0.  Zero 
warping  is  also  assumed,  although  this  is  strictly  true  only  for  beams  of 
circular  cross-section. 

Using  Timoshenko  kinematics,  the  nonlinear  strain  increments  are 
substituted  into  the  second  term  of  equation  (3.5.1),  and  after  some 
algebra,  the  volume  integral  can  be  reduced  to 
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Timoshenko  kinematics: 


5e :  aj  dV  =  j  Nx(5exoexo  +  5vc 


^  ^\\t 


+  Mz{5Kz8xo  +  5exoKz  -  5KtWo,x  -  5Wo,xKt) 

+  MyfoKyExo  +  SexoKy  +  5KtVo.x  +  5Vo,xKt) 
+  Pz(6KzKz  +  5KtKt)  +  Py(5KyKy  +  5KtKt) 
Pyz(^^y^z  SKzKyj 

+  Vy(-  50z£xo  ■  SExo^z  SGxWq^x  SWq^xQx) 

+  Vz(50y£xo  +  S£xo0y  ■  S0xVo,x  "  SVo,x0x) 

+  Ryy(5Kz0z  50zKz  "t"  SiCtOx  §0xK(j 
+  Ryz(-  5Kz0y  -  50yKz)  +  Rzy(5Ky0z  +  50zKy) 
+  Rzz(-  5Ky0y  -  60yKy  +  5Kt0x  +  60x  kJ)  dX 


(3.5.2) 


A  similar  expression  can  be  obtained  by  substituting  Bemoulli-Euler 
kinematics  into  equation  (3.5.1).  Alternatively,  the  Bemoulli-Euler  result 
can  obtained  directly  from  the  Timoshenko  result,  equation  (3.5.2),  by 
letting  Vy  =  Vz  =  0,  0z— >  Vo,x»  0y  — >  -  Wo,x> 
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Equations  (3.5.2)  &  (3.5.3)  include  all  terms.  The  higher  order 
resultants,  however,  are  easily  computed  only  for  simple  loadings.  One 
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instance  where  they  can  be  computed  is  in  the  application  of  an  axial  load. 
In  the  buckling  problem,  inclusion  of  the  Py,  Pyz,  and  P^  resultants  lowers 
the  buckling  load  for  thick  beams  only  very  slightly,  so  that  these  terms  are 
reasonably  ignored  (see  tabulations  in  Appendix  C). 

It  is  again  convenient  to  define  a  vector  of  quantities  which  will 
allow  definition  of  what  shall  be  called  the  matrix.  Let 

CL'  &  ~  [  ^*0  ^o.x  Wo,x  Ky  ] 

C^:  &  ~  [  ^*0  ^0,X  Wq  X  Kt  Ky  Kz  0X  0y  0z  ] 

The  matrix  ^  relates  the  above  quantities  to  the  element  nodal 
displacements,  and  is  formed  by  the  appropriate  placement  of  interpolation 
functions  and  their  denvatives. 

&  =  (3.5.4) 

Equation  (3.5.4)  is  valid  for  and  formulations.  Rewriting 
equation  (3.5.2)  &  (3.5.3)  in  matrix  form  and  introducing  equation  (3.5.4) 
allows  us  to  write  the  geometric  stiffness  term  in  the  following  form 


Kct  =  geometric  stiffness  matrix 

Dg  =  matrix  of  total  accumulated  stresses 

The  above  form  of  the  geometric  stiffness  matrix  is  valid  for  both 
Cl  and  CO  formulations,  with  substitution  of  appropriate  expressions  for  G 
and  Dg  is  given  below,  where  the  partition  made  up  of  the  first  six 

rows  and  coluirms  is  appropriate  to  the  Cl  formulation.  The  full  matrix  is 
the  CO  form. 

Nx  0  0  0  My  Mz  0  Vz  -Vy 

0  Nx  0  My  0  0  -  Vz  0  0 

0  0  Nx  -Mz  0  0  Vy  0  0 

OMy-MzO  0  0  0  0  0 

My  00000000 

Mz  00000000 
O-VzVyO  0  0  0  0  0 

Vz  00000000 
-Vy  00000000 

The  G.  matrix  follows  directly  from  equation  (3.5.4).  It  is  given 
explicitly  below: 


0  0  0 


0  0  0.0  0 


0  0  0  0  0  0  0  0 


0  0  n1„  0  JN^, 

0  0  ON?,  0 


0  0  0  0  JNi., 


0  0  0  0  N? 


0  0  N{,„  0  JNi,„  0  0  0  Ni,„  0  JNl, 

0  N|,j,  0  0  0  JNi.„  0  N^,„  0  0  0  JN‘ 
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For  completeness,  the  derivation  of  consistent  nodal  loads  is 
considered.  The  left  hand  side  of  equation  (3.1.1)  is  rewritten  below 


LHS  =  5q^^f£  dV  +  dS 


(3.6.1) 


70 


where,  as  a  reminder. 


U  =  flexible  displacement  of  material  particle,  with  respect  to 
the  body  fixed  frame 

fp  =  force/imit  volume,  with  respect  to  the  body  fixed  frame 
tp  =  surface  traction  applied  over  SaF»  with  respect  to  the  body 
fixed  frame 

Vp  =  reference  configuration  of  the  flexible  body 

It  should  be  noted  that  equation  (3.6.1)  is  referred  to  the  entire 
flexible  body.  With  the  introduction  of  Timoshenko  beam  kinematics,  the 
volume  integral  associated  with  the  body  force  can  be  written 

j  St|'^  If  dV  =  J  [{5uo  -  y50z  +  z60y)  fx  +  (Svq  -  z60x)  fy  +  (5wo  +  y 50x)  fz  ]  dV 

(3.6.2) 

Integration  through  the  cross-section  along  with  the  assumption  that 
beam  cross-section  principle  axes  are  aligned  with  the  beam  coordinate 
system  yields 

ri 

[  5uo  6vo  5wo  ]  fy  dx 

(3.6.3) 

It  is  easily  verified  that  Bemoulli-Euler  kinematics  produce  the  same 
equation  (3.6.3).  If  we  now  consider  the  integral  in  (3.6.3)  as  a  sum  of 
integrals  over  all  element  domains,  all  that  remains  is  to  express 
[5uo  6vo  5wo]  in  terms  of  the  element  nodal  DOFs. 


f 

J  \ 


^  dV  =  I  (5uofx  +  Svofy  +  Swofi  )  dx  = 
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where 
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Thus,  the  volume  integral  associated  with  body  forces  becomes 


Sn'^fEdV  =  5a'^ 


f 

J  \ 


Ne^fp  dx  =  5q^QB 


where 

Ne^fp  dx 

(3.6.4) 

The  surface  integral  is  dealt  with  in  the  same  manner.  After 
substitution  of  tractions  into  the  surface  integral,  an  expression  can  be 
derived  just  as  for  the  body  forces.  Equation  (3.6.1)  can  then  be  written  as 

LHS  =  89’'(%  +  Qs)  =  8aTQ  (3.6.5) 
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where 


Q  =  Qb  +  Qs 

and  comparison  with  earlier  sections  shows  that  Q  is  the  element  consistent 
nodal  load  vector. 


The  external  loads  on  a  structure  can  also  be  specified  in  an 
‘inconsistent’  fashion.  External  loads  can  be  lumped  directly  at  the  nodes. 
Direct  nodal  loading  can  be  handled  as  an  additional  term  added  to  the 
external  virtual  work.  No  interpolation  functions  are  involved  in  the 
discretization  since  the  load  is  applied  directly  to  the  element  nodal  DOFs. 
In  pr-ictice,  nodal  loads  are  added  directly  to  the  corresponding  element  of 
the  r.isembled  force  vector. 
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Chapter  4 

Eigenvalue  Problems 


A  set  of  problems  has  been  addressed  to  examine  the  beam  kinematic 
assumptions  and  finite  element  approximations.  For  a  beam  with  cantilever 
boundary  conditions,  the  following  problems  have  been  solved:  free 
vibration,  and  static  and  dynamic  buckling. 

A  study  of  the  beam  finite  elements  derived  in  the  previous  chapter 
has  been  accomplished  for  two  beams  whose  physical  and  material 
parameters  are  indicated  in  table  4.1.  The  first  beam  has  L/h  =  10 
(length/thickness)  and  borders  the  BemouUi-Euler  assumptions.  Parameter 
set  2  differs  only  in  the  beam  thickness;  this  beam  has  L/h  =  100  and  its 
behavior  should  be  consistent  with  Bemoulli-Euler  beam  theory.  This 
choice  of  parameters  allows  clearer  insight  into  the  effects  of  kinematic 
assumptions  and  treatment  of  the  mass  matrix  and  level  of  integration. 

Table  4.1.  Beam  material  properties  for  eigenvalue  problems. 


Parameter 

Parameter 

Setl 

Set  2 

E 

107 

107 

V 

0.3 

0.3 

_ P _ 

1 

1 

10 

10 

'wSSSm 

1 

1 

h 

1 

0.1 

I 

1/12 

1/12  (0.1)3 

Im 

1/12 

1/12  (0.1)3 
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Considering  only  the  flexible  domain,  the  virtual  work  expression 
gives  rise  to  the  following  dynamic  equilibrium  equation 

+  Kg  =  Q 

where 


M  =  consistent  mass  matrix 
K  =  stiffness  matrix 

Q  =  consistent  nodal  loads  (may  include  lumped  loads) 

Solutions  to  the  homogeneous  problem  can  be  obtained  by  letting 
q  =  where  y  is  the  column  matrix  of  amplitudes.  Substitution  leads 

to  a  general  eigenvalue  problem 


(-  (i)2m  +  K)j:  =  0 


(4.1.1) 


where  the  eigenvalue  is  the  square  of  the  natural  frequency  of  vibration, 
and  y  the  associated  eigenvector.  The  exact  solution  for  natural  vibration 
of  a  BemouUi-Euler  beam  is  given  by  [32] 


<0  =  (If  (Elf'" 


LI  'm 


(4.1.2) 


where  X  is  found  from  solution  of  the  characteristic  equation 


cos  X  cosh  A,  +  1  =  0 


(4.1.3) 
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The  natural  modes  of  vibration  are  shown  in  figure  4.14  and  are  given 
analytically  by 

Y(x)  =  A  [(sin  A,  -  sinh  A,)  (sin  A,x/L  -  sinh  A.x/L) 

+  (cos  A.  -  cosh  A,)  (cos  A.x/L  -  cosh  A,x/L)]  (4.1.4) 


where 

A  = _ C, _ 

sin  A,  -  sinh  A, 

Equation  (4.1.2)  yields  an  infinite  number  of  frequencies,  however 
the  higher  frequencies  obtained  become  increasingly  invalid  as  the  kinetic 
energy  associated  with  rotation  of  the  beam  cross-sections  becomes 
significant.  Recall  that  the  Bemoulli-Euler  assumptions  ignore  rotatory 
inertia;  kinetic  energy  comes  only  from  transverse  deflections.  The 
frequencies  obtained  from  equation  (4.1.2)  also  neglect  shear,  which  is 
important  as  frequency  increases,  i.e.  as  the  ratio  of  wavelength/thickness 
decreases. 

To  assess  the  effects  of  each  variation,  the  natural  frequency  was 
normalized  by  the  exact  result  from  Bemoulli-Euler  beam  theory,  and 
plotted  against  the  number  of  elements  comprising  the  beamt  (in  some 
literature  convergence  is  plotted  against  number  of  DOFs  -  equivalent). 
The  following  eight  types  of  beam  elements  have  been  formulated: 

1 )  C®:  consistent  mass,  full  integration  of  mass  and  stiffness 

2)  C^:  consistent  mass,  reduced  integration  of  stiffness 


t  A  slightly  different  perspective  is  given  in  Appendix  E,  where  the  natural 
frequencies  are  normalized  by  the  exact  result  from  Timoshenko  beam 
theory. 
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3)  C^:  consistent  mass,  1 -point  integration  of  mass  and  stiffness 

4)  C®:  lumped  mass,  full  integration  of  mass  and  stiffness 

5)  C^:  lumped  mass,  reduced  integration  of  stiffness 

6)  C^:  lumped  mass,  1 -point  integration  of  mass  and  stiffness 

7)  C^:  consistent  mass,  full  integration  of  mass  and  stiffness 

(with  rotatory  inertia  -  Rayleigh) 

8)  C^:  consistent  mass,  full  integration  of  mass  and  stiffness 

(without  rotatory  inertia  -  Bemoulli-Euler) 

Recall  that  the  consistent  mass  and  material  stiffness  matrices  were 
defined  in  equations  (3.4.9)  and  (3.4.10),  respectively.  The  consistent 
mass  and  material  stiffness  matrices  were  defined  in  equations  (3.3.15)  and 
(3.3.16),  respectively. 

By  convention,  BemouUi-Euler  theory  neglects  rotatory  inertia  and 
corresponds  directly  to  case  8.  Although  consistent  application  of  the  finite 
element  method  leads  to  fully  integrated  consistent  mass  and  stiffness 
matiices,  in  practice  a  reduced  level  of  integration  is  used  to  evaluate  the 
CO  stiffness  matrix  to  avoid  element  locking.  Reduced  integration  is 
applied  to  all  terms  of  the  stiffness  matrix,  and  is  differentiated  from 
‘selective  reduced  integration,’  in  which  only  the  shear  related  terms  are 
evaluated  using  reduced  integration.  Appendix  F  provides  a  brief 
overview  of  Gauss  quadrature,  and  gives  the  integration  rule 
corresponding  to  full  and  reduced  integration  of  the  beam  finite  elements. 
Use  of  the  lumped  mass  matrix  is  normally  dictated  by  the  reduced  cost  of 
computing  eigenvalues. 

Case  1  is  the  consistent  application  of  the  finite  element  method  to 
the  beam  element,  and  he  convergence  is  shown  in  figure  4.1. 
Generally,  natural  frequencies  converge  very  slowly,  even  for  the  first 


77 


mode.  The  increasing  influence  of  shear  and  rotatory  inertia  with  mode 
number  can  also  be  seen,  as  the  higher  modes  do  not  converge  to  the  B.E. 
result.  When  the  L/h  ratio  is  increased  by  an  order  of  magnitude,  shear 
and  rotatory  inertia  effects  are  virtually  eliminated,  and  the  problem  of 
shear  locking  becomes  more  apparent,  as  shown  in  figure  4.2. 


2D  C-zero  Cantilever  Beam:  Consistent  Mass, 
Full  Integration,  Parameter  Set  1 


#  Elements 


Figure  4.1.  Convergence  for  beam,  consistent  mass,  full  integration. 
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2D  C-zero  Cantilever  Beam:  Consistent  Mass, 
Full  Integration,  Parameter  Set  2 


>> 
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Oh 


#  Elements 


Figure  4.2.  Convergence  for  C®  beam,  consistent  mass,  full  integration. 


Case  2  is  shown  in  figure  4.3,  and  demonstrates  the  effect  reduced 
stiffness  has  on  the  C^  beam  element.  The  first  mode  is  essentially 
captured  with  one  element.  Again  the  influence  of  shear  with  increasing 
frequency  is  apparent,  but  it  can  also  be  seen  that  good  convergence  is 
achieved  with  a  relatively  coarse  mesh.  With  parameter  set  2,  the  shear 
and  rotatory  inertia  effects  are  negligible,  and  convergence  is  achieved 
with  a  fine  mesh,  as  shown  in  figure  4.4.  Locking  is  siill  observed  in  the 
higher  modes. 


79 


Normalized  Natural  Frequency 
(Bernoulli-Euler) 


2D  C-zero  Cantilever  Beam:  Consistent  Mass, 
Reduced  Stiffness,  Parameter  Set  1 


#  Elements 


Figure  4,3.  Convergence  for  beam,  consistent  mass,  reduced  stiffness. 
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2D  C-zero  Cantilever  Beam:  Consistent  Mass, 
Reduced  Stiffness,  Parameter  Set  2 


#  Elements 


Figure  4.4,  Convergence  for  beam,  consistent  mass,  reduced  stiffness. 


The  convergence  plots  associated  with  case  3  are  very  similar  to 
those  associated  with  case  2  for  the  higher  modes.  Mode  i  now  converges 
from  above  and  requires  at  least  four  elements  before  locking  on  to  the 
Bemoulli-Euler  result.  These  results  are  given  as  figures  4.5  and  4.6. 
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Normalized  Natural  Frequency 
(Bernoulli-Euler) 


Chapter  4:  Eigenvalue  Problems 


2D  C-zero  Cantilever  Beam:  Consistent  Mass, 
1 -Point  integration,  Parameter  Set  2 


#  Elements 


Figure  4.6.  Convergence  for  beam,  consistent  mass,  uniform  1 -point 

integration. 


Convergence  for  case  4  is  shown  in  figure  4.7.  The  lumped  mass 
assumption  (no  lumped  inertia)  isolates  the  shear  effect.  Comparison  with 
figure  4.1  indicates  that  rotatory  inertia  has  little  influence  on  the 
frequencies  of  vibration.  The  result  for  parameter  set  2  is  shown  in  figure 
4.8,  and  is  very  similar  to  the  consistent  mass  result  shown  in  figure  4.2, 
which  indicates  that  locking  is  not  remedied  by  mass  assumptions. 
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Normalized  Natural  Frequency 
(Bernoulli-Euler) 


2D  C>zero  Cantilever  Beam:  Lumped  Mass, 
Full  integration,  Parameter  Set  1 


#  Elements 


Figure  4,7.  Convergence  for  beam,  lumped  mass,  full  integration. 


Chapter  4:  Eisenvalue  Problems 


2D  C-zero  Cantilever  Beam:  Lumped  Mass, 
Full  Integration,  Parameter  Set  2 
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Figure  4.8.  Convergence  for  beam,  lumped  mass,  full  integration. 


The  last  set  of  C®  formulations  to  be  discussed  is  that  of  lumped 
mass/reduced  stiffness.  (Cases  5  &  6  are  one  and  the  same,  since  by  the 
lumped  mass  assumption,  the  mass  matrix  is  not  subject  to  integration). 
Note  that  only  the  effect  of  shear  is  present,  lumped  inertia  being 
neglected.  One  striking  feature  is  the  convergence  of  Mode  1  from  below, 
rather  than  from  above.  Convergence  is  achieved  with  relatively  course 
mesh,  with  shear  still  apparent  in  the  higher  modes.  The  effect  of  L/h  is 
clearly  shown  in  figure  4.10,  where  shear  is  not  an  important  factor,  and 
convergence  to  the  B.E.  result  is  observed. 
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20  C-zero  Cantilever  Beam:  Lumped  Mass, 
Reduced  Stiffness,  Parameter  Set  1 


Figure  4.9.  Convergence  for  beam,  lumped  mass,  reduced  stiffness. 


2D  C-zero  Cantilever  Beam:  Lumped  Mass, 
Reduced  Stiffness,  Parameter  Set  2 


#  Elements 


Figure  4.10.  Convergence  for  beam,  lumped  mass,  reduced  stiffness. 


The  last  two  cases  are  those  of  the  elements.  Both  are  fully 
integrated,  and  use  the  consistent  mass  matrix.  The  only  difference  is  the 
treatment  of  rotatory  inertia.  The  Rayleigh  formulation  using  parameter 
set  1  is  shown  in  figure  4. 1 1 ,  The  influence  of  rotatory  inertia  can  be  seen 
to  lower  the  frequency  of  the  higher  modes.  The  so-called  Bemoulli- 
Euler  beam  theory  has  Imy  =  Imz  =  0-  With  regard  to  the  previous 
derivation,  the  last  three  terms  on  the  diagonal  of  matrix  ^  are  set  equal 

to  zero  .  It  can  be  seen  (figure  4.12)  that  convergence  is  very  rapid,  and 
all  modes  are  converged  with  no  more  than  eight  elements.  Bemoulli- 
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Normalized  Natural  Frequency 
(Bernoulli-Euler) 


Euler  and  Rayleigh  beam  theories  become  indistinguishable  as  L/h 
increases.  Figure  4.13  is  representative  of  both  Rayleigh  and  Bernoulli- 
Euler  formulations  for  parameter  set  2. 


2D  C-one  (Rayleigh)  Cantilever  Beam: 
Consistent  Mass,  Full  Integration, 
Parameter  Set  1 


#  Elements 


Figure  4.11.  Convergence  for  beam  (Rayleigh),  consistent  mass,  full 

integration. 
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Normalized  Natural  Frequency 
(Bernoulli-Euler) 


2D  C-one  (Bern-Euler)  Cantilever  Beam: 
Consistent  Mass,  Full  Integration, 
Parameter  Set  1 


#  Elements 


Figure  4.12.  Convergence  for  beam  (B.E.),  consistent  mass,  full 

integration. 
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2D  C-one  (Rayleigh  &  Bern-Euler)  Cantilever 
Beam:  Consistent  Mass,  Full  Integration, 
Parameter  Set  2 
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Figure  4.13.  Convergence  for  beam,  consistent  mass,  full  integration. 


4.1.1  Importance  of  Shear  and  Rotatory  Inertia 

The  results  show  that  consideration  of  shear  and  rotatory  inertia 
cause  a  reduction  of  natural  frequency  compared  to  the  analytical 
Bemoulli-Euler  theory.  The  influence  of  shear  and  rotatory  inertia  can  be 
quantified,  to  a  degree,  in  terms  of  the  ‘effective  length’/thickness,  which 
allows  all  modes  and  beam  parameters  to  be  judged  according  to  the  same 
criteria.  The  effective  length  is  taken  to  be  the  approximate  wavelength, 
obtained  from  figure  4.14.  Table  4.2  shows  an  estimation  of  the 
wavelength  and  the  corresponding  ratios  of  wavelength/thickness.  For  the 
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purpose  of  discussion,  modes  1  &  2  were  considered  to  be  1/4  and  1/2 
waves,  respectively. 


Table  4.2.  Effective  (wavelength/thickness). 


Mode 

Wave¬ 

length 

Effective  length/thick 

Setl 

Set  2 

1 

40 

40 

400 

2 

20 

20 

200 

3 

8 

8 

80 

4 

5.5 

5.5 

55 

5 

4.75 

4.75 

47.5 

Correlation  of  the  values  in  table  4.2  with  the  results  given  in  the 
preceding  figures  shows  that  departure  from  BemouUi-Euler  theory  occurs 
when  the  effective  length/thickness  <  40.  This  is  demonstrated  in  figures 
4.9  &  4.10,  where  rotatory  inertia  effects  are  absent;  it  is  also  seen  in 
figure  4.11,  where  only  the  shear  effect  is  neglected.  Note  that  shear  has  a 
larger  influence  on  natural  frequency  than  does  rotatory  inertia.  As  a  rule 
of  thumb,  the  consideration  of  effective  length/thickness  provides  an  aid  in 
the  selection  of  an  appropriate  finite  element  model. 


4.1.2  Mesh  Estimate  for  CO  Elements 

An  estimate  of  the  mesh  requirements  for  CO  element  discretization 
can  be  made  with  the  help  of  figure  4.14.  The  figure  represents  the  exact 
mode  shapes  for  the  first  five  frequencies  of  transverse  vibration  for  a 
BemouUi-Euler  beam  with  cantilevered  boundary  conditions.  Amplitudes 
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shown  are  arbitrary.  Since  elements  can  model  only  linear  variation  of 
transverse  displacement,  the  number  of  straight  line  segments  (or  elements) 
required  to  model  a  given  frequency  can  be  estimated.  Discretization  based 
on  this  method  of  mesh  estimation  suggests  using  beam  elements  of  various 
lengths,  in  order  to  efficiendy  capture  the  target  mode  shapes. 


Vibration  Mode  Shapes  for  Bernouiii-Euier 
Cantiiever  Beam 


X 


Figure  4.14.  Exact  vibration  mode  shapes,  Bemoulli-Euler  theory. 
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The  effect  of  large  initial  stress  and  contribution  of  nonlinear  strains 
must  be  considered  in  order  to  solve  the  static  buckling  problem.  The 
geometric  stiffness  matrix  was  derived  in  section  3,5,  and  now  provides  the 
mechanism  by  which  the  second  eigenvalue  problem  is  formed.  The  new 
equation  can  be  written 

(K  +  Piyx  =  0  (4  2.1) 

where  the  only  nonzero  element  of  1^  is  some  reference  axial  stress 
resultant  Nxref-  The  critical  buckling  load  is  then  Per  =  A,Nxref.  where  X  is 
the  smaUest  eigenvalue.  The  exact  solution  for  BemouUi-Euler  beam  with 
clamped/free  boundary  conditions  is  given  by  [331 

PcT  =  ®^2n-lf  n  =1,2,3,... 


4L' 

where  the  associated  mode  shapes  are  given  by 

n  =  1 


(4.2.2) 


Ci(l  -  cos 

2L 


V  = 


I 


cijl  -  cos 


n  =  2,  3,  4,  ... 


(4.2.3) 


where  ci  is  arbitrary  as  long  as  the  deflections  are  consistent  with  the 
theory  of  small  displacements. 


Because  the  static  buckling  problem  does  not  involve  the  mass 
matrix,  and  because  the  load  is  constant  over  the  length  of  the  beam,  only 
three  different  formulations  are  presented  for  this  example: 


1 )  C^:  full  integration  of  stiffness  matrices 

2)  Cf*:  reduced  integration  of  material  stiffness 

3)  C^:  full  integration  of  stiffness  matrices 


Buckling  loads  have  been  normalized  by  the  exact  Bemoulli-Euler 
buckling  load.  Although  only  the  critical  mode  is  important  in  the  simple 
example  considered,  the  higher  modes  may  become  important  in  more 
complex  structures,  and  are  thus  included  in  the  following  figures. 
Convergence  of  static  buckling  load  is  very  similar  to  the  convergence  of 
natural  frequencies  obtained  from  consistent  mass  formulations  -  compare 
figures  4.15-4.20  to  figures  4.1-4. 6  and  figures  4.12-4.13. 

Comments  made  earlier  also  apply  here,  however,  a  couple  of 
additional  comments  are  made.  Figure  4.16  shows  very  clearly  the 
problem  of  shear  locking  associated  with  the  fully  integrated  C^  material 
stiffness  matrix.  Figures  4.17  and  4.18  differ  from  figures  4.3  and  4.4 
only  in  that  the  critical  buckling  is  not  captured  with  one  element.  It 
requires  at  least  eight  elements  to  lock  on  to  the  Bemoulli-Euler  result. 
The  C 1  formulation  captures  the  critical  buckling  load  with  one  element. 
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Normalized  Buckling  Load 


mhlviXi. 


Normalized  Buckling  Load 


2D  C-zero  Cantilever  Beam:  Reduced  Stiffness,  Parameter  Set  2 


#  Elements 


Figure  4.18.  Convergence  for  C®  beam,  reduced  integration  of  stiffness. 
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Normalized  Buckling  Load 


4.3  Dynamic  Buckling 


Figure  4.21.  Cantilevered  beam  under  distributed  axial  load. 


Yet  a  third  eigenvalue  problem  can  be  formulated  which  involves  all 
three  principal  finite  element  matrices.  Consider  an  inextensible  beam  with 
cantilevered  boundary  conditions  which  is  free  to  translate  in  space.  A 
constant  axial  load  is  applied  to  the  cantilevered  end  as  shown  in  figure 
4.21  which  produces  the  following  distributed  axial  loading. 


P(X)  =  p„jl  -  ij 


(4.3.1) 


For  a  given  static  loading  determined  by  Pq,  the  frequencies  of 
natural  vibration  can  be  computed.  It  will  be  seen  that  compressive  axial 
stress  has  a  ‘softening’  effect.  The  eigenvalue  problem  is  similar  to  the 
case  of  free  vibration  considered  in  section  4.1.  Now  the  effect  of  internal 
force  is  included  in  the  stiffness  and  the  eigenvalue  problem  can  be  written 
as 
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-  +  (K  +  ly]  j;  =  0 


(4.3.2) 


where  the  internal  force  distribution  is  embedded  in  the  geometric  stiffness 
matrix.  When  Pq  =  0,  the  geometric  stiffness  matrix  becomes  the  null 
matrix  and  the  simple  free  vibration  problem  is  recovered.  For  P©  <  Pcr» 
the  eigenvalue  problem  of  equation  (4.3.2)  gives  the  frequencies  of  natural 
vibration.  As  P©  is  increased  to  Per,  the  fundamental  frequency  goes  to 
zero,  as  shown  in  figure  4.22.  Dynamic  buckling  occurs  at  Pq  =  Pcr»  which 
corresponds  to  the  zero  frequency.  This  is  equivalent  to  a  static  buckling 
problem  of  the  type  considered  in  section  4.2  in  which  the  axial  load 
distribution  is  given  by  equation  (4.3.1). 

Buckling  of  a  vertical  cantilevered  beam  due  to  its  own  weight  was 
considered  by  Timoshenko  [34]  and  is  here  used  as  the  exact  solution  for 
the  Bemoulli-Euler  beam.  The  critical  load  was  given  as 

P„  =  7.837  El 

(4.3.3) 

A  more  general  formulation  for  the  buckling  of  a  beam  under  axial 
acceleration  with  rigid  mass  attached  to  the  free  end  was  considered  by 
Storch  and  Gates  [35].  In  the  degenerate  case,  with  zero  tip  mass,  the 
critical  load  has  the  same  form  as  given  above  except  the  factor  becomes 
7.8664. 

An  investigation  of  the  performance  of  a  single  element  is  shown  in 
figure  4.22.  Data  displayed  in  the  figure  is  tabulated  in  table  4.3. 
Examination  of  the  data  shows  that  the  mass  matrices  from  B.E.  and 
Rayleigh  theory  do  indeed  yield  different  frequencies,  although  the 
difference  is  very  small.  A  comparison  of  performance  would  be 
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incomplete  without  considering  the  fully  integrated  beam,  however,  a 
single  element  is  so  stiff  that  plotted  on  the  scale  of  figure  4.22  the  curve 
would  begin  well  up  (171)  and  end  well  to  the  right  (196,033)  of  all  the 
other  curves.  Reducing  the  integration  of  the  stiffness  matrices  reduces  the 
locking  problem  so  apparent  with  full  integration.  The  curve  with  reduced 
stiffness  compares  quite  favorably  to  the  Bemoulli-Euler  result. 
Overrelaxation  is  observed  in  the  C®  formulation  with  both  reduced 
stiffness  and  lumped  mass  approximation  for  small  axial  loads. 


Axial  Acceleration  of  Beam: 
Single  Element  Performance 


Acceleration 


Figure  4.22.  Natural  frequencies  for  Pq  <  Pcr»  parameter  set  1. 
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Figure  4.23  shows  the  convergence  of  the  critical  acceleration  for 
several  beam  element  formulations.  Bemoulli-Euler  and  Rayleigh  beam 
assumptions  (B.E.  has  zero  rotatory  inertia)  produce  identical  results  for 
critical  acceleration,  since  this  corresponds  to  a  static  buckling  problem 
with  load  distribution  given  by  equation  (4.3.1).  The  C®  formulation  with 
full  integration  shows  clearly  the  element  locking  problem,  while  reduced 
integration  of  the  stiffness  matrices  produces  better  convergence.  The  loss 
of  monotonicity  of  convergence  due  to  reduced  integration  is  also 
observed. 


Dynamic  Buckling:  Axial  Acceleration  Applied 
to  C-zero  and  C-one  Beams 


Figure  4.23.  Convergence  for  dynamic  buckling,  parameter  set  1. 
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Table  4.3.  Data  displayed  in  figure  4.22. 


Fundamental  Frequency  j 

Acx:cleration 

BemouUi- 

Euler 

Rayleigh 

CP,  reduced 
stiffness 

C9,  reduced 
lumped  mass 

0 

32.2 

J2l 

31.4 

25.7 

500 

- 

30.9 

- 

- 

1000 

29.7 

29.6 

29 

23.7 

2000 

26.9 

26.8 

26.3 

21.5 

3000 

23.8 

23.7 

23.2 

19 

4000 

20.2 

20.1 

19.8 

16.2 

5000 

15.8 

15.8 

15.5 

12.7 

5700 

11.8 

11.7 

- 

- 

6000 

9.53 

9.51 

9.54 

7.81 

6200 

7.69 

7.68 

- 

- 

6300 

6.59 

6.57 

- 

- 

6400 

5.25 

5.24 

- 

- 

6500 

3.42 

3.42 

4.04 

3.31 

6530 

2.64 

2.64 

- 

- 

6550 

1.95 

1.95 

- 

- 

6570 

- 

- 

2.43 

1.98 

6574.1 

0.138 

0.0555 

- 

- 

6600 

- 

- 

1.18 

0.967 

6609.3 

- 

- 

0.094 

0.0769 
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Chapter  5 

Dynamic  Problems 


The  eigenvalue  problems  of  the  previous  chapter  provide  insight 
regarding  the  behavior  of  nonlinear  beam  finite  elements  in  the  context  of 
dynamics.  The  next  step  in  the  systematic  assessment  of  finite  element 
approximations  is  a  benchmark  nonlinear  dynamic  simulation.  The 
example  chosen  was  originally  proposed  by  Kane  [36,  13]  to  demonstrate 
the  physical  behavior  of  rapidly  spiiming  systems  known  as  ‘centrifugal 
stiffening.’  The  result  reported  by  Ryan  [36]  was  obtained  using  an 
assumed  modes  formulation.  Other  researchers  who  have  also  published  a 
solution  to  the  spin-up  problem  are  Simo  &  Vu-Quoc  (nonlinear  finite 
element  method)  [15],  and  Ider  &  Amirouche  (also  assumed  modes)  [16]. 
These  published  results  provide  comparison  for  the  present  solution.  The 
focus  is  on  understanding  the  influence  of  the  Coriolis  and  centrifugal 
forcing  terms,  and  the  contribution  of  the  geometric  stiffness  matrix. 

Application  of  the  dynamics  formulation  is  also  made  to  a  space 
shuttle/remote  manipulator  arm/payload  model  which  demonstrates  a 
practical  application  of  the  theory.  A  realistic  torque  is  applied  to  the 
orbiter  and  realistic  payload  mass  properties  are  used. 
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5.1  Dynamics  Simulation 

5.1.1  Specification  of  Equations  of  Motion 

An  exact  set  of  governing  equations  was  derived  in  section  2.5.  The 
succeeding  sections  introduced  two  alternative  modelling  assumptions  with 
regards  to  the  flexible  appendage  mass  distribution:  lumped  mass  (§2.6) 
and  lumped  mass/inertia  (§2.7).  In  either  case  the  stiffness  matrix  is 
provided  by  3-D  finite  elements.  The  lumped  mass/inertia  equations  of 
section  2.7  are  chosen  for  implementation,  which  eliminates  the  necessity 
of  condensing  rotational  DOFs  from  the  stiffness  matrix.  In  the  spin-up 
problem,  all  rigid  body  DOFs  are  prescribed.  For  the  orbiter/RMS 
problem,  a  rigid  body  torque  is  prescribed  while  the  remaining  rigid  body 
DOFs  are  constrained. 

To  facilitate  discussion  of  the  results  to  follow,  equation  (2.7.15)  is 
rewritten  below  along  with  identification  of  terms  appearing  in  the  forcing 
vector  associated  with  flexible  translational  DOFs. 


Rf 

(6Nxl) 


centrifugal  Coriolis 

force  force 

^  -  mito.^  u  -  mifli^  Ti  -  2miai^  ^  +  h.o.t. 

^  +  h.o.t. 


^  -  mNCO.''  u  -  mNca.''  ^  -  2mNt!i'<  ^  +  h.o.t. 

^  +  h.o.t. 


(5.1.1) 


Note  that  the  relative  displacements  q^i  are  omitted  from  the 
centrifugal  force  term.  It  is  assumed  that  flexible  deformations  are  small 
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SO  that  vehicle  geometry  is  adequately  represented  by  the  undeformed,  or 
original  configuration.  For  the  dynamics  problems  considered,  the 
undeformed  configuration  is  the  reference  configuration.  This  implies  that 
the  rigid/flex  coupling  matrix  given  by  equation  (2.7.12)  is  constant  and 
need  be  determined  only  once  in  the  original  configuration.  Also,  the 
instantaneous  vehicle  inertia  appearing  in  the  Mrr  matrix  is  constant  and 
need  be  computed  only  once  in  the  original  configuration. 


5.1.2  Incremental  Solution 

The  equations  of  motion  (2.7.7)  are  well  suited  for  solution  by  the 
second  order  Newmark  integration  scheme.  Derivation  of  Newmark 
integration  for  linear  systems  is  outlined  in  Appendix  B,  as  well  as  an 
incremental  form  which  is  required  for  the  solution  of  nonlinear  equations. 
The  incremental  form  with  modified  Newton-Raphson  iteration  has  been 
implemented  in  the  dynamic  simulation. 

From  equation  (B.2.10),  note  that  the  incremental  solution  algorithm 
requires  calculation  of  a  nodal  force  vector  corresponding  to  the  state  of 
internal  stress.  For  the  finite  element  discretization  of  stiffness,  the 
internal  force  vector  is  given  by 
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where 


B  =  strain-displacement  matrix  encountered  in  chapter  3 
(£  =  Ba) 

a  =  vector  of  stress  resultants 

Assuming  material  linearity,  the  stress  resultants  can  be  found  using 
equation  (A.7). 

5.1.3  Computer  Implementation 

The  computer  implementation  of  the  dynamics  formulation  noted 
above  proceeds  as  shown  in  figure  5.1.  For  all  problems  the  integration 
step  size  used  was  At  =  0.01,  and  the  convergence  criteria  used  was 

II AUi  <  tolerance  (5.1.3) 

where  AU*^  is  the  vector  of  incremental  displacements  corresponding  to  the 
kth  iteration  (see  Appendix  B).  Thus  the  solution  is  converged  when  the 
norm  of  AU*'  is  less  than  some  tolerance.  The  displacement  tolerance  for 
both  dynamic  problems  was  equal  to  0.000001. 
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Simulation  Flow 

Read  control  parameters 

•  start,  stop  time;  integration  time  step  At 

Initialize  simulation 

•  read  finite  element  data  (material  stiffness) 

•  calculate  Mrf,  Mfr,  Mff  for  original  (undeformed) 
configuration 

Integration  loop 

calculate  prescribed  torques/displacements 
form  the  tangent  stiffness  matrix  Kt  based  on  Uf 

iteration  loop  (Modified  Newton-Raphson) 

•  calculate  Rf 

•  calculate  Mrr  &  Rr  from  AUTOLEV  subroutine 

•  incremental  Newmark  (5  =  1/2,  a  =  1/4) 

•  test  convergence 
—  No: 

Yes:  update  system  and  exit  iteration  loop 

data  output 
check  simulation  time 
time  <  stop  time: 

time  >  stop  time:  exit  integration  loop 
End  simulation 


Figure  5.1.  Flow  diagram  for  dynamic  simulations. 


110 


5.2.1  Problem  Description 


Figure  5.2.  Spin-up  of  cantilever  beam. 


Schematic  of  the  physical  problem  is  shown  in  figure  5.2.  A  rigid 
central  body  has  a  rigidly  attached  flexible  beam.  The  hub  radius  is  small 
compared  with  the  length  of  the  beam  and  may  be  ignored.  The  central 
body  is  constrained  against  translation,  while  hub  rotation  is  a  prescribed 
function  of  time,  equation  (5.2.1)  specifying  a  smooth  transition  from  zero 
hub  motion  to  constant  angular  speed  of  6  rad/sec.  The  prescribed  rotation 
is  also  shown  graphically  in  figure  5.3. 

jo  A  [t  -  {7.5/7c)sin  (7ct/7.5)]  rad/sec  0  <  t  <  15  sec 
\6  rad/sec  15  <  t  <  30  sec 


(5.2.1) 


Figure  5.3.  Hub  motion  prescribed  by  equation  (5.2.1). 

5.2.2  Finite  Element  Model 

The  beam  has  uniform  cross-section  and  is  discretized  using 
uniformly  sized  C®  finite  elements.  Both  material  and  geometric  stiffness 
matrices  are  evaluated  by  reduced  integration  to  avoid  shear  locking. 
Material  properties  are  chosen  to  correspond  with  previously  published 
results,  and  are  given  in  table  5.1.  (k  is  the  shear  correction  factor,  p  is 
the  per  volume  density). 
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Table  5.1.  Beam  material  properties  for  spin-up  problem. 


A  schematic  of  the  finite  element  model  in  the  undeformed 
(reference)  configuration  is  shown  in  figure  5.4.  The  body  fixed  frame  is 
coincident  with  the  rigid  central  body.  Flexible  mass  distribution  is  treated 
using  the  lumped  mass/inertia  assumption.  Note  that  the  vehicle  c.m.  in 
figure  5.4  has  only  one  nonzero  component  in  the  body  fixed  frame. 
Lumped  masses  are  equally  spaced  along  the  length  of  the  beam. 


rigid  mi  m2  m3  mN-i  mw 

central 

body 


Figure  5.4.  Schematic  of  finite  element  model,  analogous  to  figure  2.1. 
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5.2.3  Spin-UD  Results  &  Discussion 

The  study  of  the  spin-up  problem  is  accomplished  through 
consideration  of  the  following  cases: 

1)  full  solution  -  including  all  forcing  terms,  geometric 
stiffness  matrix,  incremental  integration  scheme 
(Newton-Raphson  iteration) 

2)  neglect  Coriolis  forcing  term 

3)  neglect  centrifugal  &  Coriolis  forcing  terms 

4)  neglect  geometric  stiffness  matrix 

In  all  of  the  above  cases,  the  effect  of  mesh  refinement  is  also  examined. 
The  motion  shown  in  the  following  figures  corresponds  to  the  axial  and 
transverse  displacements  of  the  beam  tip,  relative  to  the  body  fixed  frame. 
Transverse  displacement  is  perpendicular  to  the  axis  of  rotation. 

Tip  displacements  for  case  1  are  shown  in  figures  5.5-5.12.  The 
gross  behavior  is  captured  well  by  a  course  mesh,  as  seen  in  figures  5.5  and 
5.6.  Axial  displacement  appears  to  settle  into  a  steady-state  elongation  at 
constant  spin-rate.  Increasing  the  number  of  elements  brings  out  the  detail 
of  this  motion,  however,  and  it  is  seen  that  an  oscillatory  motion  is 
superposed  onto  the  steady-state  deflection.  Steady-state  axial  deflection 
converges  from  above;  in  other  words,  a  course  mesh  overpredicts  the 
axial  deflection. 

The  transverse  behavior  of  the  beam  is  characterized  by  an  increase 
in  deflection  until  the  midpoint  of  the  spin-up  (corresponding  to  the 
maximum  angular  acceleration),  after  which  the  tip  catches  up  to  the  body 
fixed  frame  and  oscillates  about  zero  relative  deflection.  This  behavior  is 
also  captured  with  a  course  mesh.  Two  tendencies  are  noted  with  respect 
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to  mesh  refinement:  the  peak  deflection  (absolute  value)  converges  from 
above,  and  the  steady  state  frequency  associated  with  bending  vibration 
converges  from  below.  For  future  reference,  case  1  is  summarized  in  table 
5.2. 
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Transverse  Displacement 


Spin-up  Problem^  2  elements 

Tip  motion 


Figure  5.5.  Axial  displacement  of  beam  tip  for  case  1,  2  elements. 


Spin-up  Problem,  2  elements 

Tip  motion 


Figure  5.6.  Transverse  displacement  of  beam  tip  for  case  1,  2  elements. 
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Transverse  Displacement 


Spin-up  Problem,  4  elements 

Tip  motion 


Figure  5.7.  Axial  displacement  of  beam  tip  for  case  1,  4  elements. 


SpIn-up  Problem,  4  elements 

Tip  motion 


Figure  5.8.  Transverse  displacement  of  beam  tip  for  case  1,  4  elements. 
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Transverse  Displacement 


Spin-up  Problem,  8  elements 

Tip  motion 


Figure  5.9.  Axial  displacement  of  beam  tip  for  case  1,  8  elements. 


Spin-up  Problem,  8  elements 

Tip  motion 


Figure  5.10.  Transverse  displacement  of  beam  tip  for  case  1,  8  elements. 
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Spin-up  Problem,  16  elements 

Tip  motion 


Figure  5.11.  Axial  displacement  of  beam  tip  for  case  1,  16  elements. 


Spin-up  Problem,  16  elements 

Tip  motion 


Figure  5.12.  Transverse  displacement  of  beam  tip  for  case  1,  16  elements. 


119 


Chapters:  Dynamic  Problems 


Table  5.2.  Summary  of  spin-up  problem,  case  1. 


#of 

elements 

Peak  transverse 
displacement 

Steady-state  axial 
displacement 

Frequency, 
bending  vibration 
(rad/sec) 

2 

-  0.753 

5.78  X  10-4 

3.36 

4 

-  0.660 

5.31  X  10-4 

3.56 

8 

-  0.646 

5.18  X  10-4 

3.77 

16 

-  0.587 

5.15  X  10-4 

3.94 

The  effect  of  neglecting  the  Coriolis  forcing  term  is  considered  in 
case  2  and  displayed  in  figures  5.13  &  5.14.  One  can  see  that  the  Coriolis 
term  provides  the  excitation  of  the  axial  vibration  mode;  its  removal 
produces  a  true  steady-state  axial  elongation.  The  Coriolis  term  is  not  a 
major  contributor  to  the  transverse  behavior  of  the  beam  (transverse 
displacements  corresponding  to  figures  5,13  &  5.14  are  the  same  as  given 
for  case  1  results). 
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Axial  Displacement  Axial  Displacement 


Spin-up  Problem,  8  elements 

Tip  motion,  no  Corioiis  term 


Figure  5.13.  Axial  displacement  of  beam  tip  for  case  2,  8  elements. 


Spin-up  Problem,  16  elements 

Tip  motion,  no  Coriolis  term 


Figure  5.14.  Axial  displacement  of  beam  tip  for  case  2,  16  elements. 
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Case  3  (no  centrifugal  or  Coriolis  forcing  terms)  results  are  shown 
in  figures  5.15-5.18.  Elimination  of  centrifugal  forcing  eliminates  any 
axial  elongation  of  the  spinning  beam.  No  other  effect  is  observed.  The 
transverse  displacement  is  unaffected  by  the  presence  or  absence  of  either 
centrifugal  or  Coriolis  terms,  as  shown  in  the  following  figures. 


Spin-up  Problem,  2  elements 

Tip  motion,  no  Coriolis  or  centrifugal 


Figure  5.15.  Transverse  displacement  of  beam  tip  for  case  3,  2  elements. 
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Time 


Omega _plot^3.0 


Figure  5.16.  Transverse  displacement  of  beam  tip  for  case  3,  4  elements. 


Spin-up  Problem,  8  elements 

Tip  motion,  no  Corioiis  or  centrifugal 


Figure  5.17.  Transverse  displacement  of  beam  tip  for  case  3,  8  elements. 
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Spin-up  Problem,  16  elements 

Tip  motion,  no  Coriolis  or  centrifugal 


Figure  5.18.  Transverse  displacement  of  beam  tip  for  case  3,  16  elements. 

Solution  of  the  spin-up  problem  in  the  absence  of  the  geometric 
stiffness  matrix  (case  4)  gives  a  surprising  result.  It  is  interesting  that  the 
geometric  stiffness  matrix  apparently  has  no  effect  on  the  solution  of  the 
spin-up  problem.  The  formulation  as  implemented  has  no  dependence  on 
geometric  nonlinearity;  material  stiffness  is  sufficient  to  prevent  divergent 
behavior.  Different  modelling  assumptions  and  solution  techniques  lead  to 
dynamic  instability.  Using  assumed  modes  formulations,  Ryan  [13]  and 
Ider  &  Amirouche  [16]  report  divergent  transverse  displacements  when 
geometric  stiffening  effects  are  ignored.  Simo  &  Vu-Quoc  [15],  on  the 
other  hand,  do  not  discuss  any  divergent  behavior  encountered  with  their 
finite  element  solution. 
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5.3  Orbiter/Remote  Manipulator  Arm  (RMS) 

5.3.1  Problem  Description 

As  a  realistic  application  of  the  dynamic  formulation,  an  additional 
example  is  considered  which  involves  the  space  shuttle  remote  manipulator 
system  (RMS).  The  manipulator  arm  is  assumed  to  be  locked  in  a  fully 
extended  configuration  with  a  typical  payload  attached  to  the  end  effector. 
The  shuttle  is  constrained  against  translational  and  rotation,  except  for 
rotation  about  the  body  fixed  y  axis.  Orbiter  rotation  is  induced  by 
application  of  positive/negative  pulse  torque,  shown  graphically  in  figure 
5.19.  The  torque  acts  about  the  y  axis  of  figure  5.20  so  that  the  elbow  is  in 
the  plane  of  rotation.  Angular  displacement  of  the  orbiter  in  response  to 
the  torque  is  heavily  influenced  by  the  relative  size  of  the  inertias  of  the 
shuttle  and  the  flexible  (appendage  +  payload)  combination. 

Motions  are  assumed  to  be  small,  thus  the  reference  configuration  is 
always  taken  to  be  the  original  configuration. 
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RMS  Problem 
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Figure  5.19.  Torque  time  history  applied  to  orbiter. 

5.3.2  Finite  Element  Model 

The  finite  element  representation  of  this  system  is  shown  in  figure 
5.20,  where  the  shuttle  c.m.  is  located  at  node  1  and  payload/end  effector 
are  located  at  node  16.  The  orbiter  and  payload  are  modelled  by  lumping 
mass  and  inertia  at  their  respective  nodes,  values  used  are  given  below. 
The  body  fixed  frame  is  coincident  with  the  orbiter  c.m.  Fifteen  elements 
are  used  to  model  the  manipulator,  which  can  be  numbered  consecutively 
starting  from  node  1.  The  finite  element  model  of  the  remote  manipulator 
arm  was  extracted  and  simplified  from  a  NASTRAN  model  (in-house 
Draper  model)  The  details  of  the  finite  element  model  are  provided  in 
Appendix  D. 
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Figure  5.20.  Finite  element  model  of  remote  manipulator  arm  in  straight 

out  position. 


morbiter  =  6397  slugs 

mpayload  —  665  slugs 
1044 

Iorbiter=  134  xlO^  slug-in^ 

1003.7  . 

10.4 

Ipayload  =  10.4  xlO^  slug'in^ 

3.469  - 


Orbiter  response  is  obtained  assuming  a  rigid  vehicle  to  provide  a 
baseline  for  the  flexible  solution,  and  provide  a  greater  sense  that  the 
flexible  solution  is  correct.  The  flexible  response  should  (and  does) 
oscillate  about  the  rigid  solution. 
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RMS  Problem 

Orbiter  motion 


ire  5.21.  Orbiter  response  to  pulse  torque,  assuming  rigid  and  flexible. 


RMS  Problem 

Orbiter  motion 


Figure  5.22.  Orbiter  response  to  pulse  torque,  assuming  rigid  and  flexible. 
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The  relative  motion  of  the  payload/end  effector,  in  the  body  fixed 
frame,  is  shown  in  figures  5.23  &  5.24,  It  is  seen  that  moderate  torque 
input  results  in  payload  deflections  on  the  order  of  1.5  inches.  This  may  be 
significant  in  the  context  of  assembling  space  station  components.  Axial 
displacements  are  not  large.  The  phase  plane  shown  in  figure  5.25 
demonstrates  the  periodic  oscillation  experienced  by  the  payload  after  the 
torques  are  released. 


RMS  Problem 

Payload  motion 


Figure  5.23,  End  effector  response  to  pulse  torque. 
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RMS  Problem 

Payload  motion 


Figure  5.24.  End  effector  response  to  pulse  torque. 


RMS  Problem 

Phase  Plane,  End-Effector 


Figure  5.25.  Demonstration  of  periodicity  of  tip  motion  (zero  damping). 
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Conclusion 


6.1  Summary  and  Conclusions 

Equations  of  motion  have  been  consistently  derived  for  a  rigid  body 
with  attached  flexible  appendage  through  application  of  the  virtual  work 
principle.  This  method  allowed  a  natural  and  consistent  introduction  of  a 
finite  element  discretization  for  the  flexible  appendage.  Using  the 
assumptions  of  lumped  mass  and  lumped  mass/inertia  allowed  comparison 
between  the  two  resulting  formulations.  The  lumped  mass  formulation 
required  condensing  the  rotational  DOFs  from  the  stiffness  matrix. 

The  nonlinear  finite  elements  used  in  conjunction  with  the  dynamic 
equations  were  consistently  derived  from  the  virtual  work  principle.  The 
consistent  derivation  using  Bemoulli-Euler  kinematics  lead  to  a  finite 
element  formulation  of  Rayleigh  beam  theory.  Bemoulli-Euler  beam 
element  was  recovered  by  setting  the  rotatory  inertias  to  zero.  All  higher 
order  terms  were  retained  in  the  derivation  of  the  geometric  stiffness 
matrix,  which  lead  to  the  introduction  of  higher  order  stress  resultants. 
Implementation  of  these  elements  allowed  greater  flexibility  than  was 
possible  with  other  commercially  available  finite  element  codes.  This 
flexibility  was  exploited  in  the  study  of  finite  element  approximations  and 
kinematic  assumptions. 

Assessment  of  element  behavior  was  accomplished  through 
eigenvalue  problems.  The  numerous  figures  generated  from  the 
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implementation  of  the  consistently  derived  finite  elements  allowed  isolation 
of  individual  effects  such  as  lumped  mass,  rotatory  inertia,  and  reduced 
integration.  Some  general  conclusions  are  drawn  from  the  results: 

•  reduced  integration  of  the  stiffness  matrix  in  elements 
alleviates  the  problem  of  shear  locking 

•  lumped  mass  assumption  lowers  vibration  frequencies 

•  shear  and  rotatory  inertia  lowers  the  vibration  frequencies 

•  shear  and  rotatory  inertia  become  increasingly  important  in 
higher  modes  of  vibration  (as  wavelength/thickness 
decreases) 

•  shear  generally  has  larger  effect  upon  frequency  compared 
to  rotatory  inertia 

•  higher  order  stress  resultants  present  in  the  geometric 
stiffness  matrix  lower  buckling  loads 

In  realistic  dynamic  analyses,  the  necessity  of  using  an  element  with  shear 
and  rotatory  inertia  was  shown  to  be  related  to  the  wavelength/thickness 
ratio.  If  the  vibration  modes  of  interest  have  large  wavelength/thickness 
ratio,  these  effects  can  be  neglected. 

The  lumped  mass/inertia  equations  of  motion  were  employed  in  the 
solution  of  two  example  problems.  In  the  spin-up  problem,  the  effects  of 
forcing  terms  and  nonlinear  flexibility  on  the  solution  were  systematically 
addressed.  It  was  shown  that  the  centrifugal  term  had  the  effect  of 
producing  a  steady-state  elongation  of  the  beam.  The  Coriolis  term  was 
responsible  for  the  axial  oscillations  superposed  on  the  steady-state 
deflection.  In  the  formulation  used  to  solve  the  spin-up  problem,  the 
material  stiffness  matrix  was  sufficient  to  prevent  divergent  behavior. 
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In  both  dynamics  problems,  the  assumption  of  small  flexible 
displacements  was  employed,  so  that  the  original  configuration  was  the 
reference  configuration.  The  formulation  as  presented  can  be  extended  to 
systems  undergoing  large  flexible  displacements  by  updating  the 
configuration.  Flexible  translations  and  rotations  are  referred  to  the 
current  configuration. 

6.2  Future  Work 

Recommended  future  work  includes  full  implementation  of  the 
rigid/flexible  appendage  formulation  to  include  configuration  updates, 
allowing  solution  of  problems  involving  large  displacements  and  rotations 
of  the  flexible  appendage.  This  involves  updating  nodal  locations,  material 
and  geometric  stiffness,  vehicle  inertia,  etc.,  at  each  time  integration  step. 
Flexible  displacements  and  rotations  are  referred  to  the  current 
configuration.  Solution  of  the  ‘spaghetti’  problem,  for  instance,  serves  as 
full  demonstration  of  the  equations  of  motion  of  chapter  2. 

The  virtual  work  principle  can  be  used  to  derive  the  governing 
equations  for  rigid  bodies  with  articulated  flexible  appendages  [28],  and 
flexible  bodies  connected  to  flexible  bodies.  These  formulations  have 
wider  application  to  spacecraft,  space  structures,  and  robotics  than  does  the 
present  formulation. 

Further  study  is  necessary  for  complete  understanding  of  the 
influence  of  geometric  nonlinearities.  Investigation  should  be  conducted 
into  the  assumed  modes  formulation  with  parallel  development  of  finite 
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element  formulation.  Systematic  study  of  this  type  should  resolve  the 
question  that  has  arisen  regarding  the  geometric  stiffness  in  the  solution  of 
dynamic  problems. 
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Figure  A.  1,  Definition  of  stress  resultants,  (a)  Stresses  at 
arbitrary  point,  (b)  Direction  of  positive  resultants. 


Stress  resultants  provide  a  convenient  measure  of  the  internal  force 
in  the  beam  elements.  Positive  stresses  are  shown  in  figure  A.l.  Stress 
resultants  are  then  defined  in  the  following  way: 
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The  strains  are  defined  as  the  work  conjugates  to  the  stress 
resultants.  The  corresponding  strains  are  therefore 
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Appendix  A:  Stress  Resultants 


QO 

C} 

^xo  ~  ^,x 

^xo  ~  Uo_x 

Ki  =  0X.X 

7\ 

II 

o 

X 

X 

K 

<x> 

1 

II 

Ky  =  Wo,xx 

~  0Z,X 

Kz  =  Vo.xx 

Yxzo  ~  '^o,x  0y 

N 

CD 

1 

K 

O 

> 

II 

o 

Let  the  vector  of  stress  resultants  be  defined  by 
Ch  a'f  =  [  N.  M,  M,  M.  ] 

CO;  =  [  N,  M,  My  Mj  V,  V,  ] 


The  stress  resultants  are  related  to  the  element  strains  by  the 
appropriate  material  matrix.  Strains  have  already  been  shown  to  be  related 
to  the  element  nodal  DOFs  by  the  B.  matrix.  The  stress  can  then  be 
evaluated  at  any  point  T|  within  the  element  from  the  equation 


=  DkB 
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Appendix  B 

Time  Integration  Schemes 


For  direct  integration  of  linear  equations  of  motion,  implicit 
integration  schemes  have  an  advantage  over  explicit  schemes  in  that  they 
are  unconditionally  stable.  The  only  restriction  on  time  step  is  due  to 
consideration  of  solution  accuracy,  and  even  then  it  need  not  be  chosen 
such  that  the  highest  frequency  is  integrated  accurately.  Newmark 
integration  is  a  second  order  implicit  scheme,  which  makes  it  an 
appropriate  choice  for  solution  of  second  order  equations  of  motion.  The 
Newmark  method  can  also  be  extended  to  the  incremental  solution  of 
nonlinear  equations  of  motion.  As  each  increment  is  a  linear  step  over  At, 
the  stability  holds  for  nonlinear  systems  as  well.  The  incremental 
Newmark  scheme  is  implemented  in  the  dynamic  simulations.  A  brief 
derivation  of  the  Newmark  method  follows. 

For  reference,  the  fourth-order  Runge-Kutta  integration  scheme  is 
also  introduced.  Because  of  its  explicit  nature,  the  time  step  must  be 
chosen  for  stability  as  well  as  accuracy,  although  for  direct  integration,  the 
stability  requirement  is  usually  strict  enough  to  assure  an  accurate  solution. 
The  Runge-Kutta  method  is  used  to  solve  systems  of  first  order  ODEs. 
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For  a  set  of  ordinary  differential  equations,  consider  the  state  of 
equilibrium  at  time  t+At: 

M  +  C  +  K  (B.l.l) 


where  M,  C,  and  K  are  the  mass,  damping,  and  stiffness  matrices,  and  R  is 
the  vector  of  external  loads.  The  three  matrices  are  constant  for  linear 
analyses.  Introduce  approximations  for  U,  and  U  at  time  t+At 


+  [(1  -  5)  ‘ti  +  5  '■'^'u]At 

(B.1.2) 

=  ‘U  +  tlAt  +  [(1  -  a)  +  a  ‘"'"^'ulAt^ 

(B.1.3) 

L\2  /  J 

where  a  and  8  are  parameters  which  govern  the  accuracy  and  stability  of 
integration.  For  5  =  1/2  and  a  =  1/4,  the  scheme  is  second-order  accurate, 
unconditionally  stable  and  equivalent  to  the  trapezoidal  mle  (also  known  as 
the  constant-average-acceleration  method). 

Equation  (B.1.3)  can  be  rearranged  for  and  substituted  into 
equation  (B.1.2).  Now  expressions  for  and  can  be  substituted 
into  equation  (B.l.l),  which  can  be  rearranged  to  give 

K‘*'“U  =  ‘*"r  (B.1.4) 

where  K  is  the  effective  stiffness  matrix  given  by 

K  =  K  +  aoM  +  aiC  (B.1.5) 
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The  integration  scheme  developed  in  the  previous  section  can  be 
viewed  as  an  incremental  solution  method  where  successive  increments  are 
At,  2At,...,  nAt;  all  increments  are  referred  to  the  original  configuration. 
This  is  possible  for  linear  equations  of  motion  since  the  coefficient  matrices 
are  constant.  In  nonlinear  problems,  the  stiffness  matrix  is  interpreted  as 
the  tangent  stiffness  matrix,  and  accounts  for  material  stiffness  as  well  as 
geometric  stiffness  and  material  nonlinear  effects.  The  tangent  stiffness  is 
configuration  dependent  and  is  denoted  at  time  t  as  The  configuration 
dependence  of  nonlinear  problems  prohibits  referring  every  increment 
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back  to  the  original  configuration,  but  an  incremental  procedure  can  be 
used  to  propagate  the  solution  from  one  known  equilibrium  state  to 
another,  providing  the  increment  is  small. 

Assume  a  state  of  equilibrium  is  known  at  time  t  which  can  be 
written  as 

m‘u  +  c'xJ  +  ‘F  =  ‘R  (B.2.1) 

where  ‘F  is  the  internal  force  vector  and  is  a  function  of  U.  With  this 
change  of  notation,  the  equilibrium  equation  (B.1.1)  can  be  written  as 

M  +  C  +  ‘^^'F  =  (B.2.2) 

Now  it  is  assumed  that  the  equilibrium  state  expressed  in  equation 
(B.2.2)  can  be  approximated  by  a  linear  increment  from  the  equilibrium 
expressed  in  equation  (B.2.1).  The  internal  force  vector  at  time  t+At  can 
thus  be  written 

=  'F  +  'Kt  AU  (B.2.3) 

where  AU  is  a  vector  of  increment  displacements.  Substitution  into 
equation  (B.2.2)  yields 

M  +  C  +  'Kt  AU  =  '-"^‘R  -  'F  (B.2.4) 

Displacements  corresponding  to  time  t+At  can  then  be  calculated  from 
t+Atu  +  AI]  (B.2.5) 

Equations  (B.2.4)  and  (B.2.5)  form  an  approximate  incremental 
solution  to  the  set  of  nonlinear  equations.  Equality  can  be  achieved 
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Appendix  B:  Time  Integration  Schemes 


through  iteration  at  time  t+At  by  using  the  full  or  modified  Newton- 
Raphson  technique.  Full  Newton-Raphson  updates  the  tangent  stiffness  at 
each  iteration,  while  the  modified  Newton-Raphson  updates  only  at  the  start 
of  each  increment.  The  full  technique  converges  more  rapidly  but  at  the 
cost  of  forming  the  tangent  stiffness  at  every  iteration.  With  the  modified 
Newton-Raphson,  the  incremental  algorithm  can  be  stated  as 


M  +  C  +  ‘Kt  AU’^  = 

(B.2.6) 

t+At^k  _  t+At^k-1 

(B.2.7) 

where 

t.itjjO  ^  tjj 

(B.2.8) 

t+AtpO  _  tp 

(B.2.9) 

Equation  (B.2.6)  is  an  incremental  form  corresponding  to  equation 
(B.1.1),  and  the  Newmark  method  of  section  B.l  can  be  applied. 
Neglecting  the  damping  matrix  and  choosing  a  and  5  corresponding  to  the 
trapezoidal  rule  gives 

K  AU*^  =  (B.2.10) 

where 

K  =  'Kt  +  -^  M 

At^  (B.2.11) 

-  M  f-^  -  *u)  -  ^  tl  -  ‘tl) 

lAt^  At  /  (B.2.12) 
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The  current  state  can  be  obtained  from 


t+Atjjk  _  t+At^k-1 

(B.2.13) 

t+At|jk  _  4  |t+At^k  tmj 

IP 

1 

I 

At^ 

At 

(B.2.14) 

^  ^  ‘11  +  ^ 

2  2 

(B.2.15) 

where  the  start-up  conditions  are  the  same  as  given  by  equations  (B.2.8) 
and  (B.2.9).  Other  schemes  have  been  developed  for  estimating  the  initial 
iteration  conditions,  rather  than  using  the  previous  equilibrium  state  as  the 
first  estimate  of  the  next  equilibrium  state  [37]. 


The  incremental  form  of  the  Newmark  method  derived  in  this 
section  for  nonlinear  problems  can  be  reduced  to  the  result  obtained  in 
section  B.l  for  linear  systems  by  noting  that  ‘Kt  ->  K  so  that  the  internal 
force  vector  becomes 


t+Atck-l 


t+AtTTk-l 


which  can  be  used  along  with  equation  (B.2.13)  to  come  up  with  equation 
(B.1.4). 


B.3  Error  Sources  in  Newmark  Integration 

For  unconditionally  stable  integration  schemes,  the  choice  of  time 
step  At  is  governed  only  by  consideration  of  accuracy.  In  structural 
dynamics  problems,  high  frequency  modes  usually  contribute  little  to  the 
response  of  the  structure.  If  the  time  step  is  chosen  to  accurately  integrate 
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only  some  subset  of  vibration  modes,  then  the  errors  contributed  by  the 
higher  modes  are  also  assumed  to  be  small.  Solution  errors  are  confined  to 
the  upper  vibration  modes  by  appropriate  choice  of  the  time  step  At. 

Since  all  modes  are  included  in  direct  integration,  solution  errors 
stem  only  from  the  use  of  too  large  an  integration  step  size.  This  is  in 
contrast  with  modal  reduction  techniques,  which  also  invite  errors  due  to 
truncation  of  the  modal  set. 

Integration  errors  are  classified  in  terms  of  period  elongation  and 
amplitude  decay  (algorithmic  damping).  The  trapezoidal  rule 
implementation  of  the  Newmark  scheme  (5  =  1/2,  a  =  1/4)  is  second-order 
accurate  and  introduces  only  period  elongation  -  no  amplitude  decay. 
Thus,  all  frequencies  contribute  to  the  structural  response.  Amplitude 
decay  can  be  introduced  in  the  Newmark  scheme  through  alternative  choice 
of  parameters  5  &  a,  although  accuracy  is  reduced  to  first-order. 

B.4  Fourth-Order  Runae  Kutta  Integration  [38] 

The  fourth-order  accurate  Runge-Kutta  is  an  explicit  multi-step 
integration  scheme  based  on  the  Euler  method.  It  operates  on  systems  of 
first  order  ordinary  differential  equations,  in  contrast  with  the  Newmark 
scheme,  which  operates  directly  on  the  second  order  equations  of  motion. 
Thus  some  additional  manipulation  is  required  in  order  to  implement  the 
Runge-Kutta  algorithm. 
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Any  system  of  ordinary  differential  equations  can  be  reduced  to  a  set 
of  N  coupled  first  order  equations.  To  preserve  the  dynamics  context,  the 
notation  of  the  previous  sections  is  used  to  express  the  general  form  as 


U,,...,  Un) 
dt 


i  =  1,...,  N 


(B.4.1) 


where  Ui  is  the  ith  component  of  the  vector  U  and  the  known  function  f^  is 
the  corresponding  derivative.  Note  that  reduction  of  the  equations  of 
motion  to  first  order  form  means  that  N  =  2(6+6*Num_nodes). 

For  convenience  only,  consider  equation  (B.4.1)  with  N=l, 
Increasing  the  system  size  involves  a  straightforward  introduction  of  a  DO 
loop  over  i=l,N.  The  Euler  method  advances  the  solution  of  a  first  order 
equation  from  t  to  t+At  by  application  of  the  formula 

t+Aty  =  ^  tyj  (B.4.2) 

where  it  should  be  noted  that  the  solution  at  time  t+At  is  based  entirely  on 
information  known  at  time  t.  Thus  no  iterations  are  necessary  in  the 
explicit  algorithm.  The  accuracy  and  stability  of  the  one  step  Euler  method 
can  be  improved  by  introducing  multiple  steps.  The  break-even  point  is 
the  fourth-order  Runge-Kutta,  which  makes  use  of  four  steps,  and  is  given 
by 


=  ‘u  +  ^^+^+^+  — 

6  3  3  6 


(B.4.3) 
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where 

ki  =  At/  '(t,  ^U) 
k2  =  At/  'jt  + 

k3  =  At/'jt  +  ^,  ‘U  +  l|j 

k4  =  At/  (t  +  At,  ‘U  +  k3) 


(B.4.4) 

(B.4.5) 

(B.4.6) 

(B.4.7) 


As  noted  earlier,  use  of  Runge-Kutta  requires  some  manipulation 
that  the  Newmark  method  does  not.  This  section  shows  the  details  of  this 
manipulation.  The  equations  of  motion  (2.6.6)  or  (2.7.7)  can  be  rewritten 
as  a  system  of  matrix  equations 

Mrr  Ur  +  Mrf  tip  =  Rr  +  Rrp  (B.5.1) 


MpR  Ur  +  Mpp  lip  +  Kpp  Up  =  Rp 


(B.5.2) 


Equation  (B.5.2)  can  be  rewritten  as 
tip  =  Mpp-'  [Rp  -  MpR  Ur  -  Kpp  Up] 


(B.5.3) 


Substitution  into  equation  (B.5.1)  yields 

[Mrr  -  Mrp  Mpp-'  Mpr]  Ur  =  Rr  +  Rrp  -  Mrp  Mpp '  Rp  +  Mrp  Mpp  '  Kpp  Up 

(B.5.4) 
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A  set  of  first  order  equations  can  be  formed  by  letting  Ur  =  Ur  and 

Up  =  Up-  The  set  of  derivatives  (at  time  t)  to  be  evaluated  in  the  course  of 

the  Runge-Kutta  scheme  can  be  summarized  as  follows; 

Ur  =  Ur  (B.5.5) 

Up  =  Up  (B.5.6) 

[Mrr  -  Mrf  Mpp'^  Mpr]  Ur  =  Rr  +  Rrf  -  Mrf  Mff  *  Rp  +  Mrf  MpF  '  Kfp  Uf 

(B.5.7) 

Mff  Up  =  Rf  -  MpR  Ur  -  KpF  Uf  (B.5.8) 

Both  of  equations  (B.5.5)  &  (B.5.7)  encompass  six  scalar  equations 
governing  the  rigid  body  DOFs.  Equations  (B.5.6)  &  (B.5.8)  encompass 
(6*Num_nodes)  scalar  equations  relating  the  flexible  DOFs.  Thus  the  total 
number  of  derivatives  evaluated  at  each  step  is  2(6  +  6*Num_nodes). 

The  Runge-Kutta  method  is  well  suited  for  nonlinear  analysis  since  it 
is  a  natural  incremental  scheme;  each  increment  is  referred  to  the  previous 
equilibrium  state.  Nonlinear  solutions  are  obtained  by  interpreting  the 
Kff  matrix  as  the  tangent  stiffness  matrix  ‘Kj,  evaluated  at  time  t.  A 

f 

typical  step  (evaluation  of  the  functions  f\  )  of  the  Runge-Kutta  algorithm 
in  the  dynamics  simulation  is  as  follows: 


151 


Given  initial  conditions  (at  time  t):  Ur,  Uf,  Ur,  Uf 

1 )  calculate  tangent  stiffness  based  on  Up 

2)  evaluate  equations  (B.5.5)  &  (B.5.6)  for  Ur,  Up 

3)  evaluate  equation  (B.5.7)  for  Ur 

4)  evaluate  equation  (B.5.8)  for  Up 

5)  use  these  results  in  calculation  of  ki’s 

A  typical  increment  implements  this  sequence  four  times  and  employs 
equation  (B.4.3).  The  tangent  stiffness  matrix  need  only  be  calculated  at 
the  beginning  of  each  increment. 
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Appendix  C 

Convergence  Data 


Convergence  data  is  tabulated  for  the  figures  shown  in  chapter  4. 


C.1  Free  Vibration 

As  a  reminder,  the  following  cases  were  considered: 

1)  C^:  consistent  mass,  full  integration  of  mass  and  stiffness 

2)  C^:  consistent  mass,  reduced  integration  of  stiffness 

3)  C9:  consistent  mass,  1 -point  integration  of  mass  and  stiffness 

4)  C^:  lumped  mass,  full  integration  of  mass  and  stiffness 

5)  C^:  lumped  mass,  reduced  integration  of  stiffiiess 

6)  CO;  lumped  mass,  1 -point  integration  of  mass  and  stiffness 

7)  Cl:  consistent  mass,  full  integration  of  mass  and  stiffness 

(with  rotatory  inertia  -  Rayleigh) 

8)  Cl:  consistent  mass,  full  integration  of  mass  and  stiffness 

(without  rotatorv  inertia  -  Bemoulli-Euler) 


Table  C.l.  Free  vibration,  case  1,  parameter  set  1. 


Mode 

Numel=1 

Numel=2 

Numel-4 

Numel=8 

Numel=16 

Numel==32 

B.  E. 

1 

171 .0 

101.0 

58.5 

40.3 

34.2 

32.5 

32.1 

2 

6820.0 

687.0 

375.0 

247.0 

207.0 

196.0 

201.1 

3 

6860.0 

1 100.0 

678.0 

553.0 

520.0 

563.2 

4 

7040.0 

2310.0 

1310.0 

1030.0 

953.0 

1 103.7 

5 

6870.0 

2130.0 

1600.0 

1470.0 

1824.5 

6 

7100.0 

31 10.0 

2260.0 

2040.0 

■■ 

7510.0 

4170.0 

3000.0 

2660.0 

8 

7980.0 

5070.0 

3800.0 

3310.0 

9 

6880.0 

4660.0 

3980.0 

1  0 

7140.0 

5570.0 

4680.0 

Table  C.2.  Free  vibration,  case  2,  parameter  set  1. 


31 .4 


5920.0 


32.6 


331.0 


4400.0 


6680.0 


32.1 


222.0 


745.0 


1820.0 


4810.0 


6100.0 


6650.0 


6900.0 


Numel=8 

NumeU16 

Numel=32 

31.9 

31.9 

31.9 

199.0 

194.0 

192.0 

560.0 

521.0 

512.0 

1110.0 

972.0 

940.0 

1870.0 

1530.0 

1450.0 

2830.0 

2170.0 

2020.0 

3940.0 

2890.0 

2630.0 

4970.0 

3680.0 

3280.0 

6850.0 

4530.0 

3950.0 

7000.0 

5440.0 

4640.0 

Table  C-3.  Free  vibration,  case  3,  parameter  set  1. 


36.3 


6830.0 


34.2 


549.0 


6930.0 


7570.0 


32.5 


245.0 


1020.0 


4190.0 


6920.0 


7300.0 


9800.0 


Numel=8 

Numel:=16 

Numel=32 

32.0 

31.9 

31.9 

203.0 

195.0 

193.0 

591.0 

527.0 

513.0 

1240.0 

994.0 

945.0 

2280.0 

1580.0 

1460.0 

4020.0 

2290.0 

2040.0 

6700.0 

3130.0 

2680.0 

6910.0 

4120.0 

3350.0 

7390.0 

5290.0 

4070.0 

7820.0 

6590.0 

4830.0 

Table  C.4.  Free  vibration,  case  4,  parameter  set  1. 


93.5 


468.0 


57.3 


328.0 


825.0 


1340.0 


40.1 


240.0 


634.0 


1150.0 


1720.0 


2270.0 


2730.0 


3030.0 


NumeU16 

NumeU32 

34.2 

32.5 

207.0 

198.0 

554.0 

530.0 

1020.0 

979.0 

1570.0 

1510.0 

2160.0 

2100.0 

2760.0 

2720.0 

3360.0 

3360.0 

3930.0 

4000.0 

4450.0 

4640.0 

B.  E. 


32.1 


201.1 


563.2 


1103.7 


1824.5 


32.1 


201.1 


563.2 


1103.7 


1824.5 


B.  E. 


32.1 


201 .1 


563.2 


1 103.7 


1824.5 
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4 


Table  C.5.  Free  vibration,  case  5,  parameter  set  1. 


Model  Numel=1 


29.9 


229.0 


31 .4 


193.0 


562.0 


1210.0 


31.8 


194.0 


524.0 


982.0 


1530.0 


21 10.0 


2630.0 


3000.0 


NumeU16 

NumeU32 

31.9 

31.9 

194.0 

194.0 

522.0 

521.0 

967.0 

965.0 

1500.0 

1490.0 

2080.0 

2080.0 

2680.0 

2700.0 

3280.0 

3340.0 

3860.0 

3970.0 

4390.0 

4610.0 

Table  C.6.  Free  vibration,  case  6,  parameter  set  1. 


29.9 


229.0 


31 .4 


193.0 


562.0 


1210.0 


31.8 


194.0 


524.0 


982.0 


1530.0 


2110.0 


2630.0 


3000.0 


Numel=16 

Numel=32 

31.9 

31 .9 

194.0 

194.0 

522.0 

521.0 

967.0 

965.0 

1500.0 

1490.0 

2080.0 

2080.0 

2680.0 

2700.0 

3280.0 

3340.0 

3860.0 

3970.0 

4390.0 

4610.0 

B.  E. 


32.1 


201.1 


563.2 


1 103.7 


1824.5 


B.  E. 


32.1 


201 .1 


563.2 


1103.7 


1824.5 


Table  C.7.  Free  vibration,  case  7,  parameter  set  1. 


Mcxte  Numel=1  Numel=2 

1  32.2  32.0 

2  310.0  200.0 

3  661.0 


Numel=4 

32.0 

199.0 

550.0 

1060.0 

1900.0 

2920.0 

4370.0 

6130.0 


Numel=8 

Numel-16 

Numel=32 

B.  E. 

32.0 

32.0 

32.0 

32.1 

198.0 

198.0 

198.0 

201.1 

546.0 

546.0 

546.0 

563.2 

1050.0 

1040.0 

1040.0 

1 103.7 

1680.0 

1670.0 

1670.0 

1824.5 

2440.0 

2410.0 

2410.0 

3300.0 

3250.0 

3240.0 

4220.0 

4150.0 

4140.0 

5630.0 

5110.0 

5090.0 

6800.0 

6120.0 

6090.0 
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Table  C.8.  Free  vibration,  case  8,  parameter  set  1. 


Model  Nurnehl  |  Numel-2  |  Numel»4  j  Numel=8  |  Numel-16 1  Numel=32 


32.2 


318.0 


32.1 


203.0 


686.0 


1990.0 


32.1 


201.0 


568.0 


1120.0 


2080.0 


3340.0 


5300.0 


8700.0 


32.1 

32.1 

32.1 

201.0 

201.0 

201.0 

564.0 

563.0 

563.0 

1110.0 

1100.0 

1100.0 

1840.0 

1830.0 

1820.0 

2760.0 

2730.0 

2730.0 

3880.0 

3810.0 

3810.0 

5160.0 

5080.0 

5070.0 

7250.0 

6540.0 

6510.0 

9200.0 

8190.0 

8140.0 

Table  C.9.  Free  vibration,  case  1,  parameter  set  2. 


53.7 


67900.0 


30.7 


215.0 


67900.0 


68000.0 


15.7 


105.0 


325.0 


661.0 


67900.0 


68000.0 


68000.0 


Numei>8 

NumeU16 

Numel=32 

7.9 

4.1 

2.2 

50.5 

25.6 

13.9 

147.0 

72.3 

39.0 

306.0 

144.0 

76.6 

541.0 

243.0 

127.0 

863.0 

373.0 

191.0 

1240.0 

538.0 

269.0 

1580.0 

741.0 

362.0 

67900.0 

989.0 

469.0 

68000.0 

1280.0 

592.0 

Table  C.IO.  Free  vibration,  case  2,  parameter  set  2. 


1.0 


58800.0 


1.0 


11.7 


39600.0 


66100.0 


28.4 


121.0 


22500.0 


51500.0 


63400.0 


67500.0 


NumeU8 

NumeU16 

1.0 

1.0 

HIEEHiil 

19.9 

18.3 

43.9 

36.8 

86.0 

63.2 

163.0 

99.2 

331.0 

147.0 

913.0 

210.0 

14400.0 

294.0 

32600.0 

406.0 

1.0 


6. 


17.9 


35.3 


58.9 


88.9 


126.0 


170.0 


222.0 


283.0 


B.E. 


32.1 


201 .1 


563.2 


1103.7 


1824.5 


B.  E. 


1.0 
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Appendix  C:  Convergence  Data 


Table  C.l  1.  Free  vibration,  case  3,  parameter  set  2, 


Mode 

Numel=1 

Numel=2 

Numel=4 

Numel=:8 

NumeU16 

Numel=32 

B.  E. 

1 

1.2 

1.1 

1.0 

1.0 

1 .0 

1 .0 

1.0 

2 

67900.0 

19.5 

8.2 

6.8 

6.5 

lllllllllllll^^ 

3 

67900.0 

41 .3 

21.1 

18.5 

18.0 

17.8 

4 

68000.0 

324.0 

49.7 

37.8 

35.5 

34.9 

5 

67900.0 

109.0 

65.9 

59.5 

57.7 

6 

68000.0 

255.0 

106.0 

90.2 

Bi 

68100.0 

773.0 

161.0 

128.0 

8 

69100.0 

4320.0 

240.0 

175.0 

9 

67900.0 

354.0 

230.0 

1  0 

68000.0 

525.0 

296.0 

Table  C.l 2.  Free  vibration,  case  4,  parameter  set  2, 


Numel=16 

Numel=32 

4.1 

2.2 

25.3 

13.8 

70.5 

38.7 

138.0 

75.7 

226.0 

125.0 

336.0 

186.0 

466.0 

259.0 

616.0 

344.0 

784.0 

441.0 

967.0 

549.0 

Table  C.l 3.  Free  vibration,  case  5,  parameter  set  2. 


Numel=l 
1 .0 
6.4 

17.8 
35.2 

58.8 
89.1 
127.0 
174.0 
232.0 
305.0 


Numel=32 

B.E. 

1.0 

1.0 

17.8 

17.8 

34.9 

34.9 

57.8 

57.7 

86.6 
121 .0 
162.0 


209.0 


262.0 


Table  C.8.  Free  vibration,  case  8,  parameter  set  1. 


32.1 


203.0 


686.0 


1990.0 


32.1 


201.0 


568.0 


1120.0 


2080.0 


Numei=8 

Numel»16 

NumeU32  | 

32.1 

32.1 

mEsmi 

201.0 

201.0 

201.0 

564.0 

563.0 

563.0 

1110.0 

1100.0 

1100.0 

1840.0 

1830.0 

1820.0 

2760.0 

2730.0 

2730.0 

3880.0 

3810.0 

3810.0 

5160.0 

5080.0 

5070.0 

7250.0 

6540.0 

6510.0 

9200.0 

8190.0 

8140.0 

Table  C.9.  Free  vibration,  case  1,  parameter  set  2. 


Mode  Numel^l  Numel-2  NumeU4  Numel>8  Numel>16  Numel=32 


1  53.7  30.7  15.7 


215.0 


67900.0 


68000.0 


105.0 


325.0 


661.0 


67900.0 


68000.0 


68000.0 


.9 


50.5 


147.0 


306.0 


541.0 


863.0 


1240.0 


68100.01  1580.0 


67900.0 


68000.0 


4.1 

2.2 

25.6 

13.9 

72.3 

39.0 

144.0 

76.6 

243.0 

127.0 

373.0 

191.0 

538.0 

269.0 

741.0 

362.0 

immm 

469.0 

1  1280.0 

592.0 

Table  C.IO.  Free  vibration,  case  2,  parameter  set  2. 


Mode  NumeUI  NumeU2  Numel>4  I  NumeU8  I  Numel-16  I  Numels32 


1.0 


58800.0 


1.0 


11.7 


39600.0 


66100.0 


28.4 


121.0 


22500.0 


51500.0 


63400.0 


1.0 


6.6 


19.9 


43.9 


86.0 


163.0 


331.0 


1.0 


6. 


18.3 


36.8 


63.2 


99.2 


147.0 


1.0 


6. 


17.9 


35.3 


58.9 


88.9 


126.0 


67500.01  913.0  210.0  170.0 


14400.0  294.0  222.0 


32600.0  406.0  283.0 


B.  E. 


32.1 


201.1 


563.2 


1103.7 


1824.5 


B.  E. 


1.0 


B.  E. 


1.0 


6. 


17.8 


34.9 


57.7 


3 


Table  C.l  1.  Free  vibration,  case  3,  parameter  set  2. 


Numel=4 

Numel=8 

Numel=16 

Numel=32 

1.0 

1.0 

1 .0 

1 .0 

8.2 

6.8 

6.5 

41.3 

21.1 

18.5 

18.0 

324.0 

49.7 

37.8 

35.5 

67900.0 

109.0 

65.9 

59.5 

68000.0 

255.0 

106.0 

90.2 

68100.0 

773.0 

161.0 

128.0 

69100.0 

4320.0 

240.0 

175.0 

67900.0 

354.0 

230.0 

68000.0 

525.0 

296.0 

Table  C.l 2.  Free  vibration,  case  4,  parameter  set  2. 


28.3 


146.0 


15.3 


90.1 


238.0 


415.0 


48.5 


134.0 


257.0 


417.0 


604.0 


793.0 


939.0 


B.  E. 


1 .0 


6. 


17.8 


34.9 


57.7 


Numel=16 

Numel»32 

4.1 

2.2 

25.3 

13.8 

70.5 

38.7 

138.0 

75.7 

226.0 

125.0 

336.0 

186.0 

466.0 

259.0 

616.0 

344.0 

784.0 

441.0 

967.0 

549.0 

Table  C.13.  Free  vibration,  case  5,  parameter  set  2. 
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Table  C.14.  Free  vibration,  case  6,  parameter  set  2. 


1.0 


6.3 


18.0 


36.7 


65.6 


114.0 


212.0 


559.0 


Numel=16 

Numel=:32 

1.0 

1 .0 

IfllllEQH 

HHOH 

17.8 

17.8 

35.2 

34.9 

58.8 

57.8 

89.1 

86.6 

127.0 

121.0 

174.0 

162.0 

232.0 

209.0 

305.0 

262.0 

Table  C.15.  Free  vibration,  case  7,  parameter  set  2. 


Mode  Numelal 


NumeU 


1.0 


NumeU4 1  Numel=8  |  Numel»16  |  Numel-32 


1.0 


6. 


17.9 


35.4 


65.9 


106.0 


168.0 


275.0 


17.8 

17.8 

17.8 

35.0 

34.9 

34.9 

58.0 

57.7 

57.7 

87.2 

86.3 

86.2 

123.0 

121.0 

120.0 

163.0 

161.0 

160.0 

229.0 

207.0 

206.0 

291.0 

259.0 

257.0 

Table  C.16.  Free  vibration,  case  8,  parameter  set  2. 


17.9 


35.4 


65.9 


106.0 


168.0 


275.0 


Numel=8 

Numela16 

1.0 

1.0 

6.4 

WBomi 

17.8 

17.8 

35.0 

34.9 

58.0 

57.7 

87.2 

86.3 

123.0 

121.0 

163.0 

161.0 

229.0 

207.0 

291.0 

259.0 

1.0 


6. 


17.8 


34.9 


57.7 


86.2 


120.0 


160.0 


206.0 


257.0 


B.  E. 


1.0 


Appendix  C:  Convergence  Data _ 

C.2  Static  Buckling 

Static  buckling  results  were  obtained  for  the  following  cases: 

1)  C®:  full  integration  of  stiffr'iss  matrices 

2)  C^:  reduced  integration  of  material  stiffness 

3)  C^:  full  integration  of  stiffness  matrices 

Included  among  the  tabulated  data  are  the  results  using  the  higher  order 
stress  resultants.  Table  entries  have  units  of  stress. 


Table  C.17.  Static  buckling,  case  1,  parameter  set  1. 


Mode 

Nume!=1 

Numel»2 

Numel=4 

Numel>8 

NumeU16 

Numel=32 

B.  E. 

1 

9.80E+05 

2.28E+05 

7.05E+04 

3.28E+04 

2.35E+04 

2.12E+04 

2.06E+04 

2 

2.62E+06 

6.69E+05 

2.91E+05 

2.04E+05 

1.84E+05 

1.85E+05 

3 

1.98E+06 

7.78E+05 

5.31E+05 

4.73E+05 

5.14E+05 

Bi 

3.55E+06 

1,44E+06 

9.47E+05 

8.35E+05 

1.01E+06 

[5 

2.19E+06 

1.40E+06 

1.22E+06 

1.67E+06 

Table  C.18.  Static  buckling,  case  2,  parameter  set  1. 


Mode 

Numel=1 

Numel=2 

Numel=4 

Numel»8 

Numel=16 

Numel=32 

B.  E. 

1 

3.30E+04 

2.27E+04 

2.10E+04 

2.06E+04 

2.05E+04 

2.05E+04 

2.06E  +  04 

2 

6.46E+05 

2.24E  +  05 

1.87E+05 

1.79E+05 

1.77E+05 

1.85E+05 

3 

9.11E+05 

5.26E+05 

4.70E+05 

4.57E+05 

5.14E+05 

4 

2.99E+06 

1.05E+06 

8.51  E+05 

8.11  E+05 

1 .01E+06 

5 

1.74E+06 

1.28E+06 

1.19E+06 

1.67E+06 

Table  C.19.  Static  buckling,  case  3,  parameter  set  1. 
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Table  C.20.  Static  buckling,  case  1,  parameter  set  2. 


Mode 

Numel=l 

Numel=2 

Num6l=4 

Numel=8 

Numel=16 

Numel=32 

B.  E. 

1 

9.62E+05 

2.O8E+05 

5.03E+04 

1.26E+04 

3.30E+03 

9.78E+02 

2.06E+02 

2 

2.54E  +  06 

5.00E+05 

1.16E+05 

2.99E+04 

8.81  E+03 

1.85E+03 

3 

1.65E+06 

3.40E+05 

8.39E+04 

2.45E+04 

5.14E+03 

mm 

3.44E+06 

7.15E+05 

1.67E+05 

4.83  E+04 

1.01E4  04 

1.29E+06 

2.83E+05 

8.02E+04 

1.67E+04 

Table  C.21.  Static  buckling,  case  2,  parameter  set  2. 


Mode 

Numel=1 

Numel=2 

Numel»4 

Numel»8 

Numel=16 

Numel=32 

B.  E. 

1 

3.33E+02 

2.29E+02 

2.1 1E+02 

2.07E+02 

2.06E+02 

2.06E+02 

2.06E  +  02 

2 

7.76E+03 

2.38E+03 

1.96E+03 

1.88E+03 

1.86E+03 

1.85E+03 

3 

1.19E+04 

6.09E+03 

5.35E+03 

5.19E+03 

5.14E+03 

mm 

1.30E+05 

1.43E+04 

1.09E+04 

1 .02E+04 

1 .01E+04 

5 

3.14E+04 

1.90E+04 

1.71  E+04 

1.67E+04 

Table  C,22.  Static  buckling,  case  3,  parameter  set  2. 


Mode 

NumeUi 

Numel=:2 

Numel=4 

Numel»8 

Numel»16 

Numels32 

B.  E. 

1 

2.07E+02 

2.06E+02 

2.06E+02 

2.06E+02 

2.06E+02 

2.06E+02 

2.06E+02 

2 

2.68E+03 

1.91  E+03 

1.86E+03 

1.85E+03 

1.85E+03 

1.85E+03 

1.8SE+03 

3 

6.42E+03 

5.23E+03 

5.15E+03 

5.14E+03 

5.14E+03 

5.14E+03 

BM 

1.66E+04 

1.06E+04 

1.01E+04 

1.01E+04 

1.01  E+04 

1.01E+04 

5 

1.97E+04 

1.68E+04 

1.67E+04 

1.67E+04 

1.67E  +  04 

Table  C.23,  Static  buckling,  case  1,  parameter  set  1,  higher  order  terms. 


Mode 

Numel=1 

Numel=2 

Numel=4 

Numel^S 

Numel»16 

Numel=32 

B.  E. 

1 

'9"78e7o5 

2.28E+05 

7.04E+04 

3.28E+04 

2.35E+04 

2.12E+04 

2.06E  +  04 

2 

2.60E+06 

6.59E+05 

2.86E+05 

2.01  E+05 

1 .80E+05 

1.85E+05 

3 

1.92E+06 

7.49E+05 

5.10E+05 

4.54E+05 

5.14E+05 

Bi 

3.52E+06 

1.36E+06 

8.91  E+05 

7.85E+05 

1.01E+06 

5 

2.06E+06 

1.29E+06 

1.13E+06 

1.67E+06 
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Table  C.24.  Static  buckling,  case  2,  parameter  set  1,  higher  order  terms. 


Mode 

Numel=1 

Numel=2 

Numel=4 

Numel=8 

Numel=1 6 

Numel=32 

B.  E. 

“i 

3.29E  +  04 

2.27E+04 

2.09E+04 

"2.05E+04 

2.04E+04 

2.04E+04 

2.06E+04 

2 

6.13E+05 

2.20E  +  05 

1.84E+05 

1.76E+05 

1.74E+05 

1.85E+05 

3 

8.51E+05 

5.03E+05 

4.51  E+05 

4.40E+05 

5.14E+05 

mm 

2.76E  +  06 

9.70E+05 

7.97E+05 

7.61  E+05 

1.01E+06 

5 

1.58E+06 

1.17E+06 

1.10E+06 

1.67E+06 

Table  C.25.  Static  buckling,  case  3,  parameter  set  1,  higher  order  terms. 


Mode 

Numel=1 

Numel=2 

Numel=4 

Numel=8 

Numel=16 

Numel=32 

B.  E. 

1 

2.07E+04 

2.05E  +  04 

2.05E+04 

2.05E+04 

2.05E+04 

2.05E+04 

2.06E+04 

2 

2.61E  +  05 

1.88E  +  05 

1.82E+05 

1.82E+05 

1.82E+05 

1.82E+b5 

1.85E+05 

3 

6.03E  +  05 

4.97E+05 

4.89E+05 

4.89E+05 

4.89E+05 

5.14E+05 

1.42E  +  06 

9.60E+05 

9.19E+05 

9.16E+05 

9.15E+05 

1.01E+06 

LlJ 

1.65E+06 

1.44E+06 

1.43E+06 

1.43E+06 

1 .67E  +  06 

Table  C.26.  Static  buckling,  case  1,  parameter  set  2,  higher  order  terms 


Mode 

Numel=1 

Numel=>2 

Numel=4 

Numel»8 

NumeU16 

Numel=32 

B.  E. 

1 

9.62E+05 

2.08E+05 

5.03E+04 

1.26E+04 

3.30E+03 

9.78E+02 

2.06E+02 

2 

2.54E+06 

5.00E+05 

1.16E+05 

2.98E+04 

8.81  E+03 

1 .85E+03 

3 

1.65E+06 

3.39E+05 

8.39E+04 

2.45E+04 

5.14E  +  03 

KM 

3.44E+06 

7.14E+05 

1.67E+05 

4.82E+04 

1 .01E  +  04 

LjJ 

1.29E+06 

2.83E+05 

8.00E+04 

1 .67E+04 

Table  C.27.  Static  buckling,  case  2,  parameter  set  2,  higher  order  terms. 


Mode 

Numel=1 

Numel=2 

Numel=4 

Numel=8 

Numel=16 

Numel=32 

B.  E. 

1 

3.33E+02 

2.29E+02 

2.1 1E+02 

2^07E+02 

2.06E+02 

2.06E+02 

2.06E+02 

2 

7.75E+03 

2.38E+03 

1.96E+03 

1.88E+03 

1.86E+03 

1.85E  +  03 

3 

1.19E  +  04 

6.08E+03 

5.34E+03 

5.18E  +  03 

5.14E+03 

WjfJM 

1.29E+05 

1.43E+04 

1 .09E+04 

1 .02E+04 

1.01E+04 

5 

3.13E+04 

1.90E  +  04 

1.71E+04 

1.67E+04 

Table  C.28.  Static  buckling,  case  3,  parameter  set  2,  higher  order  terms. 


Mode 

Numel=1 

Numel=2 

Numel=4 

Numel=8 

Numel=16 

Numel=32 

B.  E. 

“i 

2.07E+02 

2.06E+02 

2.06E+02 

2.06E  +  02 

2.O6E+02' 

2.06E+02 

2.06E+02 

2 

2.68E+03 

1.91E+03 

1.85E  +  03 

1.85E+03 

1.85E+03 

1.85E+03 

1.85E+03 

3 

6.42E+03 

5.23E+03 

5.14E+03 

5.14E+03 

5.14E+03 

5.14E+03 

4 

1.65E+04 

1.06E+04 

1.01  E+04 

1.01  E+04 

1.01  E+04 

1.01E+04 

5 

1.97E+04 

1.68E+04 

1.66E+04 

1.66E+04 

1.67E+04 

C.3  Dynamic  Buckling 

Results  tabulated  below  correspond  to  the  convergence  plot  shown  in 
figure  4.23.  The  exact  solution  is  that  of  Timoshenko  [34]. 


Table  C.29.  Dynamic  buckling,  parameter  set  1. 


#  Elements 

Critical  Acceleration 

C1  (B.E.) 

Ravieiqh 

0°  Full 

Exact  B.E. 

1 

6574.1 

6574.1 

196033 

6609.3 

6530.8 

2 

6547.4 

6547.4 

68564 

6957.6 

4 

6532.4 

6532.4 

22159 

6635.2 

8 

6531 .1 

6531.1 

10388 

6520.0 

1  6 

6531 .1 

6531.1 

7455 

6489.3 

32 

6531 .0 

6531 .0 

6722 

6481  .5 
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Appendix  D 

RMS  Finite  Element  Model 


Details  are  provided  for  the  simplified  finite  element  model  of  the 
remote  manipulator  arm  extracted  from  the  NASTRAN  model. 
Information  given  here  was  used  in  the  finite  element  formulations  derived 
in  chapter  3  and  subsequently  used  in  the  dynamic  formulation  of  chapter 
2.  An  example  problem  is  given  in  chapter  5. 


Table  D.l.  Nodal  data  for  RMS  model. 


Node  # 

X  (in) 

7  (in) 

z  (in) 

1 

'“""6."'"" 

0. 

0. 

2 

22.5396 

0. 

0. 

3 

34.5396 

0. 

0. 

4 

84.7495 

0. 

0. 

5 

134.9594 

0. 

0. 

6 

185.1694 

0. 

0. 

7 

235.3793 

0. 

0. 

8 

285.6156 

0. 

-6. 

9 

341.2056 

0. 

0. 

10 

396.7956 

0. 

0. 

11 

452.3856 

0. 

0. 

12 

507.9756 

0. 

0. 

13 

563.5956 

0. 

0. 

14 

581.5956 

0. 

0. 

15 

611.5512 

0. 

0. 

16 

637.5962 

0. 

0. 
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Table  D.2.  Element  data  for  RMS  model 


Mat’lSet# 


1 


1 


2 


2 


2 


2 


2 


2 


2 


2 


2 


2 


1 


1 


1 


Table  D.3.  Cross  section  Properties  for  RMS  model 


Prop. 

Set# 

Imy 

(slg-in2) 

Imz 

(slg*in2) 

(in4) 

HR 

A 

(in^) 

m 

(slg/in) 

m 

1 

.03346 

.02792 

54.91 

45.81 

26.27 

.01601 

.53 

2 

.03653 

.06155 

71.81 

121.00 

30.04 

.01528 

.53 

3 

.007259 

.007259 

66.21 

66.21 

28.84 

.003162 

.53 

— 

.003367 

.003367 

44.27 

44.27 

23.59 

.001794 

.53 

5 

.002035 

.002035 

7.22 

7.22 

9.53 

.002685 

.53 

6 

.01025 

.01025 

17.20 

17.20 

14.70 

.008762 

.53 

7  .02741 

.02741 

9.7 

9.7 

3.51 

.009917 

.53 

Table  D.4.  Material  Properties  for  RMS  model 


Mat’lSet#  E 


R&iStail 


1 

1  Xl07 

.3 

2 

2.22.x  107 

.495 

Appendix  E:  Solution  of  Timoshenko  Frequence  Equation 


Appendix  E 

Solution  of  Timoshenko  Frequency  Equation 


The  analytical  formulation  of  Timoshenko  beam  theory  leads  to 
much  more  complicated  differential  equations  than  does  classical  Bemoulli- 
Euler  beam  theory.  Kruszewski  [39]  gives  the  frequency  equation  for  a 
uniform  beam  with  cantilevered  boundary  conditions  as 

2  -  ^  sin  keP  sinh  k^a 

V 1  -  ks^  kRi^  kfi^ 


ke^  (ks^  -  kRi^f  +  2  cos  keP  cosh  kea  =  0 


(E.l) 
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along  with  the  definitions: 


m  =  mass  of  the  beam  per  unit  length 
As  =  effective  shear-carrying  area 
At  =  effective  total  cross-sectional  area 
ks  =  coefficient  of  shear  rigidity 
kRi  =  coefficient  of  rotary  inertia 
ke  =  frequency  coefficient 
0)  =  natural  frequency  (rad/sec) 

For  the  beam  used  in  the  eigenvalue  problems  of  chapter  4, 
As  =  At.  The  roots  of  the  frequency  equation  (E.l)  must  be  determined 
by  numerical  methods,  and  natural  frequencies  subsequently  found  using 
equation  (E.4).  For  the  property  sets  defined  in  chapter  4,  the  natural 
frequencies  corresponding  to  the  first  5  vibration  modes  are  obtained  from 
Timoshenko  beam  theory  and  are  shown  in  table  E.l. 


Table  E.l.  Natural  frequencies  from  Timoshenko  beam  theory. 


Mode 

Natural  Frequency 

Parameter 

Setl 

Parameter 

1 

31.88 

1.01 

2 

192.07 

6.36 

3 

508.58 

17.79 

4 

929.37 

34.83 

5 

1424.45 

57.50 

These  values  may  be  compared  with  the  Bemoulli-Euler  results 
tabulated  in  Appendix  C.  Note  that  shear  and  rotatory  inertia  have  the 
effect  of  reducing  the  natural  frequencies.  For  parameter  set  1,  the 
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frequency  of  mode  5  is  reduced  by  20%.  Almost  no  difference  is  observed 
for  parameter  set  2. 

The  natural  frequency  data  for  the  element  (parameter  set  1  only) 
is  normalized  by  the  Timoshenko  solution  and  shown  in  figures  E.1-E.5 
below.  The  analogous  figure  from  chapter  4  is  indicated,  and  serves  as 
comparison.  The  magnitude  of  the  shear  and  rotatory  inertia  effects  was 
demonstrated  earlier.  Normalization  by  the  Timoshenko  solution  validates 
that  the  finite  elements  are  producing  the  correct  result. 

The  slow  convergence  of  all  modes  using  full  integration  is  again 
noted  in  figure  E.l,  This  demonstrates  the  problem  of  ‘shear  locking.’ 
Figure  E.2  shows  the  improvement  in  mode  1  convergence  achieved 
through  the  use  of  reduced  integration  of  the  stiffness  matrix. 

Note  figures  E.4  and  E.5.  Rotatory  inertia  is  not  present  in  the 
lumped  mass  assumption.  The  convergence  to  higher  frequencies  than 
predicted  by  the  Timoshenko  solution  are  directly  attributable  to  the 
exclusion  of  rotatory  inertia  in  the  finite  element  model. 

The  Rayleigh  and  Bemoulli-Euler  finite  element  results  are  shown  in 
figures  E.6  and  E.7.  From  these  figures  it  is  seen  that,  (a)  the  frequencies 
are  higher  than  predicted  by  Timoshenko  beam  theory,  and  (b)  the 
frequencies  are  reduced  by  the  inclusion  of  rotatory  inertia  (figure  E.6). 
The  difference  remaining  in  figure  E.6  is  attributable  to  the  exclusion  of 
shear  in  the  finite  element  model. 
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Figure  E.l.  Convergence  plot  normalized  by  Timoshenko  frequencies; 

analogous  to  figure  4.1. 
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Figure  E.2.  Convergence  plot  normalized  by  Timoshenko  frequencies; 

analogous  to  figure  4.3. 


Normalized  Natural  Frequency 
(Timoshenko) 


20  C-zero  Cantilever  Beam:  Consistent  Mass, 
1 -Point  Integration,  Parameter  Set  1 


0  10  20  30  40 


#  Elements 


Figure  E.3.  Convergence  plot  normalized  by  Timoshenko  frequencies; 

analogous  to  figure  4.5. 
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Figure  E.4.  Convergence  plot  normalized  by  Timoshenko  frequencies; 

analogous  to  figure  4.7. 
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Figure  E.5.  Convergence  plot  normalized  by  Timoshenko  frequencies; 

analogous  to  figure  4.9. 
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Figure  E.6.  Convergence  plot  normalized  by  Timoshenko  frequencies; 

analogous  to  figure  4.1 1. 


» 


173 


Normalized  Natural  Frequency 
(Timoshenko) 


2D  C-one  (Bern-Euler)  Cantilever  Beam: 
Consistent  Mass,  Full  Integration, 
Parameter  Set  1 


#  Elements 


Figure  E.7.  Convergence  plot  normalized  by  Timoshenko  frequencies; 

analogous  to  figure  4.12. 
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Appendix  F 

Gauss  Quadrature  [23] 


i 


The  spatial  integrals  arising  in  the  finite  element  method  are 
conveniently  evaluated  using  numerical  integration  procedures.  An 
efficient  scheme  is  of  considerable  importance  in  reduction  of  both 
computational  time  and  cost.  Gauss  quadrature  is  such  an  optimal 
procedure.  Integrals  such  as  given  by  equations  (3.3.15),  (3.3.16),  and 
(3.5.5)  are  transformed  such  that 


/{x)dx 


(F.l) 


where  /(x)  is  representative  of  an  arbitrary  element  of  the  matrix  triple 
product  and  g(^)  is  the  transformed  integrand.  The  transformation 
involves  the  Jacobian,  which  was  discussed  in  chapter  3.  In  one  dimension, 
the  quadrature  formula  is  given  by 


1=  =  WigKi)+ W2g{y+--  + Wng(y 


(F.2) 


where 


g(^)  =  transformed  function  to  be  integrated 
=  sampling  points 
wi  =  weights 
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Appendix  F:  Gauss  Quadrature 


In  Gauss  quadrature,  the  sampling  points  and  weights  are  prescribed 
such  that  maximum  accuracy  is  achieved  for  a  given  n.  Sampling  points 
are  symmetric  about  the  center  of  the  interval  and  are  unevenly  spaced. 
Symmetric  points  have  the  same  weight.  Exact  integration  is  often 
referred  to  as  ‘full  integration’,  whereas  ‘reduced  integration’  usually 
refers  to  an  order  of  integration  one  less  than  required  for  full  integration. 
The  order  of  integration  is  governed  by  the  rule:  a  polynomial  of  degree 
(2n  - 1)  is  integrated  exactly  by  n-point  Gauss  quadrature. 

The  mass  and  material  and  geometric  stiffness  matrices,  derived  in 
chapter  3,  are  evaluated  by  Gauss  quadrature.  The  integration  order 
required  for  exact  integration  can  be  determined  from  examination  of  the 
matrix  elements  which  constitute  the  integrand,  along  with  the  above  rule. 
The  integration  rule  for  the  beam  elements  is  summarized  in  table  F.  1 . 


Table  F.l.  Integration  rule  for  beam  finite  elements. 


Element 

Mass  Matrix 

Stiffness  Matrix 

Formulation 

Full 

Full 

Reduced 

Cl 

4-point 

CO 

2-point 

1 -point 

tin  practice,  reduced  integration  is  not  necessary  in  elements. 


